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ABSTRACT 



The formation and large-scale propagation of Poynting-dominated jets produced 
by accreting, rapidly rotating black hole systems are studied by numerically integrating 
the general relativistic magnetohydrodynamic equations of motion to follow the self- 
consistent interaction between accretion disks and black holes. This study extends 
previous similar work by studying jets till t w 10 4 GM/c 3 out tors 10 4 GM/c 2 , by 
which the jet is super- fast magnetosonic and moves at a lab-frame bulk Lorentz factor 
of r ~ 10 with a maximum terminal Lorentz factor of Too < 10 3 . The radial structure 
of the Poynting-dominated jet is piece-wise self-similar, and fits to flow quantities 
along the field line are provided. Beyond the Alfven surface at r ~ 10-100GM/c 2 , 
the jet becomes marginally unstable to (at least) current-driven instabilities. Such 
instabilities drive shocks in the jet that limit the efficiency of magnetic acceleration 
and collimation. These instabilities also induce jet substructure with 3 < T < 15. 
The jet is shown to only marginally satisfy the necessary and sufficient conditions 
for kink instability, so this may explain how astrophysical jets can extend to large 
distances without completely disrupting. At large distance, the jet angular structure 
is Gaussian-like (or uniform within the core with sharp exponential wings) with a 
half-opening angle of w 5° and there is an extended component out to sa 27°. Unlike 
in some hydrodynamic simulations, the environment is found to play a negligible role 
in jet structure, acceleration, and collimation as long as the ambient pressure of the 
surrounding medium is small compared to the magnetic pressure in the jet. 

Key words: accretion disks, black hole physics, galaxies: jets, gamma rays: bursts, 
X-rays : bursts 



1 INTRODUCTION 

Gamma-ray bursts (GRBs), active galactic nuclei (AGN), 
and black hole x-ray binary systems exhibit relativistic jets 
that are believed to be driven by accreting black holes (see , 
e.g.. lRees et~ai]|l982l : iBeeelman et alJll984t IWooslevlll993) . 
Accreting black holes are the highly efficient engines neces- 
sary to account for the observed kinematics and energetics 
of such relativistic jets. Understanding systems powered by 
black hole accretion requires a general relativistic magneto- 
hydrodynamics (GRMHD) model of the bulk flow dynamics 
near the black hole where the relativistic jet is launched. 
The ideal MHD approximation has been shown to be a rea- 
sonably valid model to describe the nonradiative dynami- 
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cally important black hole accretion physics associated with 
GRB s, AGN, and black hole x-r ay binary systems (see, e.g., 
|PhinnevllT98a: lMcKinnevll2004l) . This approximation is the 
foundation of most studies of jets and winds. 

We are concerned with explaining how relativistic jets 
are formed, accelerate, collimate, and what is the resulting 
structure at large radii. The primarily difficulty has been 
to explain how observed jets come to be so well-collimated 
with opening angles of a few degrees, fast with bulk Lorentz 
factors of up to T ~ 10 3 , and why magnetized jets would not 
be disrupted by instabilities, such as the kink instability. It 
has not even been known whether the MHD approximation 
was sufficient to understand the self-consistent production of 
jets or their observed fast speeds and small opening angles. 

GRMHD numerical models of black hole accretion sys- 
tems have significantly progressed our understanding of rel- 
ativistic jets. Such simulations have displayed primarily two 
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types of "jets," where one is associated with the disk that 
is mass-loaded by disk mater ial and the other is associated 
direc tly with the black hole l| De Villier^ Hawley, fc Krolikl 

li\L 



GRMHD numerical models of winds from the inner 
part of the accretion di sk have found that the Lorentz fac- 
tor only reaches V < 3 <De Villiers. Hawlev. fc Krolikll2003l 



M< -Khun -v V ( : ; u; i i ni jl2004l: iDe Villiers et aLll2005ali ". This 
includes GRMHD numerical models of GRBs that showed a 
disk wind with a velocity of only v ~ 0.3c ijMizuno et alJ 
|2004) . Thus, disk winds may not explain jets with high 
Lorentz factors unless additional physics, such as possibly 
efficient neutrino-annihilation for GRBs or shock accelera- 
tion of particles in general, is included. 

The most astrophysically plausible mechanism to gener- 
ate highly relativistic jets is via the electromagnetic extrac- 
tion of black hole spin energy by the Bl andford-Znajek (BZ) 
mechanism ijBlandford fc ZnaiekHl977f) . They modelled the 
disk as a current sheet that sources the magnetic field that 
couples to a rotating black hole, which spins down due to a 
transfer of rotational energy of the space-time to the mag- 
netic field. However, their force- free model does not include 
matter, and so cannot directly explain the Lorentz factor of 
observed jets. Also, the BZ model does not self-consistently 
determine the current in the disk or the field geometry that 
threads the black hole. All these features require the inclu- 
sion of matter, or the inertial effects of matter, that would 
be present in an accretion disk. 

GRMHD numerical models of black hole accretion 
have shown polar regions with Poynting-dominated jets 
that are consistent with being powered by the BZ effect 
( McKinn ey fc Gamrniel 120041 : iKomissarovl 120051 : IMcKinnevI 
2005a). They have also shown that the magnetic field corre- 
sponding to the Blandford-Z najek effect domina t es all other 
magnetic field geometries jHirose et alJ 12004 IMcKinnevI 
2005a). In contrast to the disk wind that is self-consistently 
mass-loaded by the disk itself, the Poynting-dominated jet 
that is launched by the BZ effect has an arbitrary Lorentz 
factor that is determined by the detailed high-energy physics 
of the mass-loading of the jet, so lo w mass-loading can ex- 
plain highly relativistic jets (see, e.g. |p hinney||l983|: [Beskinl 
Il997l: IPunslvlEoOll iLevinsorjEool lMcKinneJ2005alfb|) . 

We study a generic black hole accretion system in the 
GRMHD approximation and consider its application to all 
black hole accretion systems. This unification scheme works 
because the GRMHD equations of motion scale arbitrarily 
with the mass of the black hole and the mass accretion rate. 
Hence, a GRMHD model can be applied to any black hole 
accretion system as a unified model, apart from important 
high-energy effects that break this scale invariance. 

The primary modifications to the GRMHD model come 
from radiative and other high-energy effects. As in any other 
study, we assume that the results of a GRMHD model are 
interesting as long as such physics plays a simple role in 
the bulk jet dynamics. This role includes setting the mass- 
loading of the Poynting jet, where in our GRMHD model 
this is set to some fixed value. Also, the jet speed may be 
slowly modified, such as by Compton scattering or by syn- 
chrotron emission of the internal energy. In this case the 
radiative effects can be treated as small corrections to the 

background state set by the GRMH D model. 

This paper extends the work of iMcKinnev fc Gamrniel 



(2004) by studying the self-consistently generated Poynting- 
dominated jet as it propagates to large dis tances. A new 
version of the GRMHD code called HARM iGammie et alJ 
2003a) is used to integrate the equations of motion. Unlike 
prior work, the jet is shown to reach bulk Lorentz factors of 
r ~ 10 and maximum terminal Lorentz factors with Too < 
10 3 and collimate to narrow half-opening angles of 0j ~ 
5°. The simulated jet is sufficiently fast and collimated to 
account for jet observations associated with GRBs, AGN, or 
x-ray binary systems. 



Outline of Paper 

In Section we discuss the astrophysical systems to which 
this study can be applied. We discuss the direct or indirect 
observations of jet speed, collimation, and structure, which 
are the principle issues this paper attempts to address. 

In Section [3] we summarize the equations solved and 
the units used in the paper. Examples are given for how to 
apply our results to any astrophysical system. 

In Section |1J we present the setup of the numerical 
model, such as the initial and boundary conditions. The as- 
trophysical applicability of the model setup to GRBs, AGN, 
and x-ray binaries is considered. 

In Section |S] we present the results of a fiducial numer- 
ical model and describe the jet formation, structure, and 
propagation to large distances. Both the radial and angu- 
lar structures of the jet are studied and data fits to flow 
quantities are given. The jet characteristic structure, such 
as the formation of a double transonic (trans-fast) flow, is 
discussed. 

In Section HJ we discuss the results and their implica- 
tions. The results are compared to similar investigations, 
and the limitations of the models presented are discussed. 

In Section |7] we provide a short summary of the key 
points. 



2 ASTROPHYSICAL MOTIVATION 

Observations associated with GRBs, AGN, and black hole 
x-ray binary systems all show direct evidence of relativistic 
outflows. While AGN and x-ray binaries show direct evi- 
dence of jet collimation, outflows associated with GRBs are 
believed to be collimated based on achromatic breaks in af- 
terglow light curves. However, the formation, structure, ac- 
celeration, and collimation of relativistic jets are not yet un- 
derstood. The following discussion outlines the motivation 
for using a GRMHD model to study jets from GRBs, AGN, 
and x-ray binary systems. 

2.1 GRBs 

Neutron stars and black holes are associated with the most 
violent of post-Big Bang events: supernovae and some GRBs 
an d maybe some x-ray flashes (XRFs) (for a general r eview 
see lWooslevll993llWheeler. Yi. Hoflich. fc Wan J200(il) . Ob- 
servations of a supernova light curve (SN2003dh) in the af- 
terglow of GRB 030329 suggest that at least some long- 
duration GRBs are consistent with core- c ollapse events 
JStanek et alJ 120031: Ikawabata et al.l l2003t lUemura et all 
120031 iHiorth et al\\20o£ ~ 
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Neutrino processes and magnetic fields are both 
important to understand core- collapse. In unraveling the 
mechanism by which core-collapse supernovae explode, 
the details of the neutrino transport model was at some 
point realized to be critic al to whether a supernova is 
prod uced in simulations (Mcss er et al. and collaborators! 
1998). This has thus far been interpreted to imply that 
highly accurate neutrino transport physics is required, 
but this could also mean additional physics, such as 
a magnetic field, could play a significant role. Indeed, 
all core-collapse events may be power ed by MHD pro- 



19701: ISvmbalistvl 


11984 
Il992l 


IWooslev & Weaver! 


Duncan & Thompson! 


iKhokhlov et alJ 


Akivama. Wheeler. Meier. & Lichtenstadtl 



1986 
1999 
2003 

Core-collapse involves shearing subject to the 
Balbus-Hawley instability as in a cc retion disks 
J Akivama. Wheeler. Meier, fc Lichtenstadtl l2003f) . and 
this process creates a dynamically important magnetic field 
during collapse. All core-collapse explosions are signifi- 
cantly polarised, asymmetric, and often bi-polar indicating 



Wang: 


fc Wheelerl 


119961: 


IWheeler. Yi. Hoflich. fc Wang 


2000: IWang. Howell. Hoflich 
2002:: IWang. Baade. Hoflich 


& Wheeler 
& WheeleiJ 


l200ll: IWane et alJ 
I2OO3L and refer- 



ences therein). Possible evidence for a magnetic do minated 
outflow has been found in GRB 021206 llCoburn fc Boggs 



2003), marginally consisten t with a magnetic outflow di- 
rectly from the inner engine jLvutikov. Pariev. fc Blandfordl 
2003), although these observations remain controversial. 
However, most core-c ollapse supernovae may still not 
require magnetic fields IWooslev et alJ BOOS). 

The key energy source for many GRB models is a black 
hole accretion system since the efficiency is high for convert- 
ing gravitational binding energy to some type of emission. 
Collapsar type models are based upon the formation of a 
black hole during the core-collapse of rapidly rotating mas- 
sive stars. The typical radius of the accretio n disk may de- 
termine the duration of long-duration G RBs (IWooslevl l993: 
|Paczvnskll998l : lMacFadven fc Wooslevll999l) . An accretion 
disk is also formed as a result of a neut ron star or black hole 
collis ions with another stellar object iNaravan et al.lll992l 
2001). 



GRBs are believed to be the result of an ultrarelativis- 
tic jet. Indirect observational evidence of relativistic mo- 
tion is suggested by afterglow achromatic light breaks and 
the "compactness problem" suggests GRB material must be 
ultrarelativistic with Lorentz factor r > 100 to e mit the 
observed nonthermal 7-rays (see, e.g., |Piranll2005l) . Direct 
observational evidence for rela tivistic motion c omes from 
radio scintillation of the ISM iGoodmanl 11997^) and mea- 
surements of the afterg low emitting region from GRB030329 
jTavlor et alJl2004alliy . 

Typical GRB jet models invoke either a hot neutrino- 
driven jet or a cold Poynting flux-dominated jet. While 
both mechanisms extract compar able amounts of the ac- 
cretion energy to power the jet dPopham et al.lll999ll . a 
neutrino-driven jet derives its energy from the annihi- 
lation of neutrinos produced by the release of gravita- 
tional binding energy, and consequently the jet is ther- 
mally accelerated. However, strong outflows can be magneti- 



cally driven feisn ovatvi-Kogan fc Ruzmaik inll976l : lLovelacel 
1976HBlandfordll976HBlandford fc Znaieklll977l) . 



2.2 AGN 

AGN have long been believed t o be powered by accretion 
onto supermassive black holes (IZel'dovichl Il964l : ISahaeterl 
Il964h . Observations of MCG 6-30-15 show an iron line fea- 
ture consi stent with emission from a relativistic disk with 
v/c ~ 0.2 llTanaka et al]ll995llFabian et al]l2002l) . although 
the lack of a temporal correlation between the continuum 
emission and iron-line emission may suggest the basic fluo- 
rescence model is incorr ect. Alterna tively, the iron-line may 
be a jet-related feature (Elvis 2000). 

AGN are observed to have jets with F ~ 10 



JUrrv & Padovanil Il995l iBiretta et alJ 


19991). maybe even 


r ~ 30 jBeeelman et al.1 119944 


Ghisellini & Celottil 



d et al. 2UU1J, wnnc some observations imply 
r < 200 iGhisellini et alJ Il993l: iKrawczvnski et alJ 120021 : 
iKonopelko et alJl2003l)~ Some radio-quiet AGN show evi- 
dence of weak iets dGhisellini et al]l2004l) . which could be 
explained as a disk wind and not require a rapidly rotat- 
ing black hole jMcKinnev fc Gammiei l2004h . Observations 
imply the existence of a two-component jet structure with 
a Po ynting jet core and a diss i pative surrounding co mpo- 
nent iGhisellini fc Madaijll996l : bhisellini et al.ll2005l) . The 
energy structure of the jet and wind are important in under- 
standing the feedback effect that controls si ze of the black 
hole and may determine the M — a relation dSpringel et alJ 
l2004UDi Matteo et al.ll2005l) . 



2.3 X-ray Binaries 

Long after their formation, neutron stars and black 
holes often continue to produce outflows and jets 
jMirabel fc Rodriguez! Il999l) These include x-ray binaries 
(for a review see lLewiriet^llT995tlMcClintock fc Remhtoxd l 
|2003j), neutron star as pul sars (for a review seelLp rj^^a 
|20(H| on ms pulsars and iThorsett fc Chakrabartvl " ^99l 
on radio pulsars) and soft-gamma ray repeaters (S GRs) 
llThompson fc Duncaji99ill996UKouveliotou et al.ll999l) . 
In the case of x-ray binaries, the companion star's stellar 
wind or Roche-lobe forms an accretion disk. Many x-ray bi- 
naries in their hard/low state (and radio-loud AGN) show 
a corre lation between the x -ray luminosity and radio lumi- 
nosity jMerloni et alJ i2003). which is consistent with radio 
synchrotron emission from a jet and x-ray emission from a 
geometrically thick, optically thin, Comptoni zing disk- 
S o me blac k hole x -ray binaries have jets iMirabel et alJ 
119921: iFended I2003al1 . such as GRS 1915+105 with 
appa r ently superluminal motion (|Mjrabel^^^todrigue2j 
Il994t iMirabel fc Rodriguez! Il999t iFender fc Bellonil 12004 



Kaiser et al J l2004h . However, it is difficult to constrain the 
Lorentz factor, which in some ejection events may be as 
large for x-ray binaries as for AGN with a Lorent z factor of 
T < 100 jFenderll2003bl: iMiller- Jones et alJl2006l) . The ob- 
served broad, shifted, and asymmetric iron line from GRS 
1915+105 is possible evi dence for a relativistic accretion disk 
iMartocchia et all2002Tl . although this feature could be pro- 
duc ed by a jet component as suggested for AGN. 

iGierliriski fc Done! <l2004f) suggest that at least some sys- 
tems, such as GRS 1915+105, have slowly rotating black 
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holes (but see lCui et aljll998h . If this is correct, then the 
BZ mechanism may not be responsible for these jets. A 
baryon-loaded disk wind with T < 3 can be produced from 
a black hole ac cretion disk and not require a rapidly rotat- 
ing black hole iMcKinnev fc Gammid I2004T) . Nonrelativis- 
tic outflo ws were found even in viscous hydrodynamic sim- 
ulations llStone e t al. Ill999t llgu menshchev fc Abramowicz I 
ll999tl200fllMcKinnev fc Gammidl2002l) . 



3 EVOLUTION EQUATIONS AND UNITS 

The model is that of a rotating black hole, described by 
the Kerr metric, surrounded by an accretion disk. The Kerr 
metric is written in Kerr-Schild coordinates, such that the 
inner-radial computational boundary can be placed inside 
the horizon and so out of causal contact with the flow. The 
Kerr metric in Kerr-Schild coordinates and the Jacobian 
transformation to Boyer-Lin dquist coordinates are given in 
iMcKinnev fc Gammid <2004) . 

Boyer-Lindquist coordinates are not chosen because it is 
difficult to avoid interactions between the inner-radial com- 
putational boundary and the jet. The coordinate singular- 
ity at the event horizon in Boyer-Lindquist can be avoided 
by placing the inner-radial computational boundary outside 
the horizon. However, Poynting-dominated flows have waves 
that propagate outward even arbitrarily close to the event 
horizon. Using Boyer-Lindquist coordinates can lead to ex- 
cessive variability in the jet since the ingoing superfast tran- 
sition is not on the computational grid, and then the details 
of the boundary condition can significantly influence the jet. 
Numerical models of viscous flows h ave historically had re- 
lated issues (see discussion in, e.g., IMcKinnev fc Gammid 
2002). 



3.1 GRMHD Equations of Motion 

The GRMHD notation follows iMisner et all Jl973ft . here- 
after MTW. A single-component MHD approximation is as- 
sumed such that particle number is conserved, 



0. 



(1) 



where po is the rest-mass density and is the 4-velocity. 
A 4-velocity with a spatial drift is introduced that is unique 
by always being related to a physical observer for any space- 
time and has well-behaved spatially interpolated values, 
which is useful for numerical schemes. This 4-velocity is 



777 



(2) 



where 7 = —u a r\ a . This additional term represents the spa- 
tial drift of the zero angular momentum (ZAMO) frame de- 
fined to have a 4-velocity of 7? M — {—a, 0,0,0}, where a = 
l/\/ — 9 U an( i so u* — 7/a. One can show that 7 = (1+g 2 ) 1 ' 2 
with q 2 = gijffv? . 

For a magnetized plasma, the conservation of energy- 
momentum equations are 



MA + J EMJ 



0. 



(3) 



where T^ v is the stress-energy tensor, which can be split 
into a matter (MA) and electromagnetic (EM) part. In the 
fluid approximation 



T^ A = {p +u g )u»u» +p g P^, (4) 

with a relativistic ideal gas pressure p g = (7 — l)w S) where 
u g is the internal energy density and the projection tensor 
is P M " = g^" + u^u" , which projects any 4- vector into the 



comoving frame (i.e. P" 



0). 



In terms of the Faraday (or electromagnetic field) tensor 



EM — r r 



1 UU riO 

I 9 



(5) 



which is written in Heaviside-Lorentz units such that a fac- 
tor of 47r is absorbed into the definition of F^", where the 
Gaussian unit value of the magnetic field is obtained by mul- 
tiplying the Heaviside-Lorentz value by \/47r. The induction 

* flu 

equation is given by the space components of F ]v = 0, 
where F^ = ■|e M " KA F KA is the dual of the Faraday ten- 
sor (Maxwell tensor), and the time component gives the 
no-monopoles constraint. Here e is the Levi-Civita tensor, 
where ^ uXS = — / 1 [fiu XS] and [pi^AiS] is the completely 
antisymmetric symbol. The comoving electric field is defined 
as 

1 n.vkX *t-i -it / f\ 

m , (6) 



UuF» 



LlukX *77i 



where 77 corresponds to a scalar resistivity for a comoving 
current density j M = J„P" M . The comoving magnetic field 
is defined as 

1 



b" 



.F 



u v Fy, 



(7) 



The ideal MHD approximation, 77 = e M = 0, is assumed, and 
so the invariant e M fe M = 0. Since the Lorentz acceleration on 
a particle is / ; M = qe M , then this implies that the Lorentz 
force vanishes on a particle in the ideal MHD approximation. 
Since e"u u = b v u v = 0, they each have only 3 independent 
components. One can show that 



b = bru 



and 



pf 



(8) 



(9) 



so that the electromagnetic part of the stress-energy tensor 
can be written as 

-(«Y + P"V!)V. (10) 



EM 



The other Maxwell equations, J M = F n " - v , define the cur- 
rent density, J M , but are not needed in the ideal MHD ap- 
proximation for the evolution of the matter or the magnetic 
field. 

For numerical simplicity, another set of field vectors are 

* it 

introduced, such that B l = F and Ei = Fu/y/—g. The 
two 4- vectors and 6 M and the 3-vectors B 1 and Ei are just 
different ways of writing the independent components of the 
Faraday or Maxwell tensors. Equation Q implies b' = B'm 
and b 1 — (B 1 + u l b t ~)j-a t . Then the no-monopoles constraint 
becomes 



(V=gB l ), t = 0, 

and the magnetic induction equation becomes 

(^TgB l ), t = -(V=ff(6V - bV)),,- 

= -(^(SV-BV), 



(11) 



(12) 
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where v' = u'/u*, Ei = £i = —e i j k v :i B k — — v x B is the 
electromotive force (EMF), and e 11 is the spatial permuta- 
tion tensor. The above set of equations are those that are 
solved. A more complete discuss i on of the relativistic MHD 
equations can be found in lAnild (I1989T) . 

3.2 Floor Model 

In GRMHD numerical models of Poynting jets, the 
black hole's polar region is alway s completely evac- 
uated jMcKinnev fc Gammid 120041: iKomissarovl |2004 
| 2005l: iDe Villiers. Hawlev. fc Krolikl 120031: iDe Villiers et all 

l2005al). Unlike in t he radi al direction through the 

disk JSpruit fc Taaml 1199(1 llgumenshchev et alJ l2003t 
iNaravan et al]l2003ft . the jet is interchange stable across 
the interface b etween the disk wind a nd the Poynting- 
dominated jet jLevinson fc Eichlerlll993T) . so these nonax- 
isymmetric modes do not lead to significant mass-loading of 
the Poynting jet. Anomalous resistivity due to reconnection 
between the wind and funnel is not strong enough to signifi- 
cantly fill the funnel even at modest resolutions that under- 
resolve the funnel-w all interface iMcKinnev fc Gammid 
120041: lMcKinnevll2004|) . Also, the disk wind, which straddles 
the Poynting jet, has a sufficiently large angular momentum 
to a void significantly contam inating the Poynting jet (see, 
e.g., iMochkovitch et aljll993T) . That is, unlike in the radial 
direction in the disk, there is no angular momentum trans- 
port across the 9 direction between the Poynting-dominated 
and the disk wind. Therefore, to evolve the GRMHD equa- 
tions of motion, which require matter, a numerical scheme 
or physical mechanism must be invoked to fill the funnel 
with matter once it evacuates. 

As done by all others, in this paper the rest-mass and in- 
ternal energy density are limited to a small but finite value in 
the Poynting-dominated jet by employing a so-called "floor 
model." This model forces a minimum on the rest-mass den- 
sity (pfi) and internal energy density (w/i), which are usu- 
ally set to several orders of magnitude lower than the disk 
density. This effectively modifies the equations of motion. 

Floor-models necessarily violate the ideal MHD ap- 
proximation because they violate rest-mass and energy- 
momentum conservation. This occurs, at least, where the 
poloidal velocity u p = at the stagnation surface in the 
Poynting-dominated jet. Matter inside the stagnation sur- 
face falls into the black hole, while matter outside it is 
ejected as part of the jet. Arbitrary floor-models violate the 
ideal-MHD condition far away from the black hole where 
no rest -mass should be added for any physica l mass-loading 
model (|Phinnevll983l:lMcKinnevl2004l2005bl) . For example, 
if a numerical model uses pfi ~ Const, and the floor value is 
reached at small radii, then the ideal-MHD approximation is 
violated for the entire length of the jet in the funnel because 
the density would otherwise have decreased with radius. 

We use a floor model with a minimum allowed rest-mass 
density of pfi = 10 _7 r _2 ' 7 po,<«sfc an d minimum allowed in- 
ternal energy density of Ufi = W~ 9 r~ 27 po,disk- The coeffi- 
cients are not chosen arbitrarily small or large for two rea- 
sons. First, the value of fe 2 /po in the jet determines the mag- 
netic energy per unit particle, which determines the maxi- 
mum Lorentz factor the jet can obtain. As long as b 2 /po > 1 
near the poles, then we find that a relativistic jet emerges 
around the poles. The important physics in the emergence of 



a BZ-driven jet is that the gradient in the toroidal magnetic 
field dominates the gas+ram pressure of the surrounding 
material. The density floor is chosen to allow the emergence 
of a fast jet since the magnetic energy per unit particle is 
large (b 2 /po < 10 7 since b 2 ~ po,diskC 2 is expected near 
the poles of the black hole) in the jet once it is evacuated. 
Conversely, if a density floor coefficient is chosen such that 
Pfi ^ 10" 1 po,disk, then we find that no Poynting-dominated 
jet emerges since the value of b 2 /po < 1 around the black 
hole near the poles. A self-consistent model of the mass- 
loading must at least determine the time-averaged density 
near the black hole around the poles. Second, the floor coef- 
ficient is chosen to be not too low in order to avoid numerical 
difficulties in integrating the equations when the floor is ac- 
tivated in regions of large magnetic energy density per unit 
rest-mass density. This is a common problem with all present 
GRMHD numerical methods since the equations of motion 
become a stiff set of equations in the limit b 2 /po 2> 1. 

The radial scaling of the floor model is chosen to avoid 
contamination in the jet by the floor model at large radii. 
At large radii the rest-mass density drops off oc r~ 2 ' 2 , so the 
rest-mass density scaling of oc r~ 2,7 does not interfere with 
this. A radial scaling with a power < 2.2 would lead to non- 
physical results by loading the jet with extra high-velocity 
material. The internal energy density floor is never reached 
in the jet due to the shock heating, which keeps u/ po > I in 
the jet. See section|K|for a more detailed discussion. 

Our code keeps track of how often the density goes be- 
low the floors and how this modifies the conservation of 
mass, energy, and angular momentum. The floor model con- 
tribution is negligible except in the Poynting-dominated jet 
within r < 10r g . This injection of mass is a simple ad- 
hoc model for the mass-loading of the Poynting jet. Beyond 
r ~ 10r 9 , the floor is not activated because the accumulated 
rest-mass is larger than the floor value. Thus, beyond this 
point the jet is self-consistently evolved strictly within the 
ideal MHD approximation. 

No explicit reconnection model is included. However, 
our code checks the effective resistivity by measuring the 
rest-mass flux across field lines. For the boundary between 
the disk wind and the Poynting-dominated jet, the total 
rest-mass flux across the field line is negligible compared 
to the rest-mass flux along the field line determined by the 
floor model. An unresolved model would load field lines with 
rest-mass that crosses field lines and not properly represent 
any physical model of resistivity jMcKinnevll2004h . 

This "floor model" prescription was found to give quan- 
titatively sim ilar results to phy sically-based models of the 
mass-loading iMcKinney|2005cri for which most of the mass- 
loading t akes place only nea r the black hole for r < 10r 9 . 
Note that IMcKinnevI l|2005ch performed identical (including 
random seed for initial condition perturbations) simulations 
as in this paper except the method used to add mass and 
energy to the Poynting jet. In this paper, the floor model 
is chosen to most closely matc h the physically-based model 
described in lMcKinnevi (|2005crl . As mentioned there, the re- 
sults are quantitatively similar except within r < 10r g and 
only inside the Poynting jet. 
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3.3 Lorentz Factors 

This paper studies the speed of the jet as determined by the 
Lorentz factor. The Lorentz factor of the jet can be measured 
with respect to any frame. We measure the Lorentz factor 
with respect to a static observer at infinity, 



t t i 

u = it v—gtt, 



(13) 



in Boyer-Lindquist coordinates, where no static observers 
exist inside the ergosphere. This is in contrast to the def- 
inition of W = u*\/— l/g u , which is the Lorentz factor as 
measured by a ZAMO observer as used by many numerical 
relativists. 

Strict upper limits on the Lorentz factor and the angu- 
lar velocity at large distances can be obtained from the local 
flow quantities. The Lorentz factor and (^-velocity at large 
distances for an axisymmetric, stationary flow have upper 
limits of 
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= hus + §B<t, 



(14) 
(15) 



where p indicates either poloidal direction (e.g., r or 9), 
h = (po + u g + p) I po is the specific enthalpy, $ = B p /(pou p ) 
is the conserved magnetic flux per unit rest-mass flux, Qf = 
Ftp/Fptj, is the conserved field rotation frequency, = F$t 

is the covariant toroidal magnetic field, and uto = is 
the orthonormal angular velocity at large distances. For a 
stationary, axisymmetric flow, these are the energy and an- 
gular momentum flux per unit rest-mass flux that are con- 
served along flux surfaces. The matter and electromagnetic 
pieces are separable, such that 



= r 



(MA) 



+ r 



(EM) 



(16) 

= C '+C- (17) 

See lMcKinnevh2005bl) for more details. 

Excluding radiative processes and with some physical 
accounting of the mass- loading of the jet, the results of this 
paper suggest one can approximately assume almost all the 
magnetic and thermal energy go into particle kinetic energy. 
Thus Too and u^> i00 may represent the Lorentz factor and 
angular velocity at large distances before the energy is lost 
to radiation once the flow is optically thin. 



3.4 Units 

The units in this paper have GM = c = 1, which sets the 
scale of length (r g = GM/c 2 ) and time (t g = GM/c s ). The 
mass scale is determined by setting the observed (model- 
dependent measured or inferred for GRB-type systems) 
mass accretion rate (Mo) equal to the accretion rate through 
the black hole horizon as measured in a simulation. So 
the mass scale is scaled by the mass accretion rate (Mo) 
at the horizon (r — th = r g (l + \/l — j 2 )), such that 
Po,disk = Mo[r = ru]t g /r g and the mass scale is then just 
m = po.disfc^g = Mo[r = ru]t g - For a black hole with angu- 
lar momentum J = jGM 2 /c, j = a/M is the dimensionless 
Kerr parameter with 

The results of the simulations can be applied to 



any astrophysical system once the value of po,disk is es- 
timated. For example, a collapsar model with M = 
O.IMqs" 1 and M w 3Mp, th en p 0tdisk « 3.4 x 10 10 gcm~ 3 
jMacFadven fe Wooslevl [l999) . M87 has a mass accretion 
rate of M ~ 10 ~ 2 M^ yr - 1 and a black hole m ass of 
M » 3 x 1O 9 M jHclll999l : [Reynolds et al.lll996ft giving 
po.disk ~ 10 _16 gcm -3 . GRS 1915+1 05 has a mass accre- 
tion r ate of M ~ 7 x lCr^Mnyr " 1 iMirabel fc Rodriguez! 
119941: iMirabel fc Rodriguez! Il999l iFender fc Bellonil 120041) 
with a mass of M ~ 14M^ (iGreiner et alJl200ll) . but see 
iKaiser et all J2004h . This gives p ,di S k ~3x 10" 4 gcm -3 . 



4 NUMERICAL SETUP 

The GRMHD equations of motion are integrated numeri- 
cally usi ng a modified version of a numerical code called 
HARM iGammie et al.ll2003al) . which uses a conservative, 
shock-capturing scheme. Compared to the original HARM, 
the inversion of conserved quantities to primitive vari- 
ables is performed by using a new faster and more r obust 
two-dimensional non-linear solver jNoble et alJ 120051). A 
parabolic interpolation scheme JColella fc Woodwardlll984l) 
is used rather than a linear interpolation scheme, and an 
optimal total variational di minishing (TVD) third order 
Runge-Kutta time stepping (|Shuj|l997) is used rather than 
the mid-point method. For the problems under considera- 
tion, the parabolic interpolation and third order time step- 
ping method reduce the truncation error significantly, in- 
cluding regions where b 2 /po 2> 1. 



4.1 Computational Domain 

The computational domain is axisymmetric, with a grid 
that extends from Vi n — 0.98rn, where th is the hori- 
zon, to r out = 10 4 r 9 , and from 8 — to = 7T. Our 
numerical integrations are carried out on a uniform grid 
with coordinates: xo, xi, X2, £3, where xo = t[Kerr — Schild], 
2:3 = <^)[Kerr — Schild]. The radial coordinate is chosen to be 



r — Ro + e x 



(18) 



where Ro is chosen to concentrate the grid zones toward 
the event horizon (as Ro is increased from to th) and 
n r is controls the enhancement of inner to outer-radial re- 
gions. For studies where the disk and jet interaction is of 
primary interest, _Ro = and n r = 1 are chosen. For this 
paper, for which in addition the far-field jet is of interest, 
Ro = —3 and n r = 10 are chosen in order to reach large 
distances with few grid points and to have the inner-radial 
computational region resolve d similarly to the models of 
iMcKinnev fc Gammid J2004:). The 6 coordinate is chosen 
to be 



9 = HX2 + — (1 — h(r)) sin(27rj:2), 



(19) 



where h(r) is used to concentrate grid zones toward the 
equator (as h is decreased from 1 to 0) or pole (as h is 
increased from 1 to 2). The jet at large radii is resolved 
together with the disk at small radii using 



h(r) = 2 



i(r/roj) 



(20) 



with the parameters of Qj = 1.3, roj = 2.8, rij = 0.3, r\j 
20, r 2j = 80, and 
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9j = 9j(r) = - + -atan( ) 

Z it rij 



(21) 



is used to control the transition between inner and outer- 
radial regions. The values of Qj, r<y, rij, rij, and T2j were 
set empirically to approximately follow the simulated jet 
from the horizon to large distances. An alternative to this 
fixed refi nement of the jet and disk is an adaptive refinement 
(see, e.g., IZhang fc MacFadverjEooih . 

The numerical resolution of all the models described 
is 512 x 256 compared to th e fiducial model with 456 2 in 
iMcKinnev fc Gammid ll2004) . However, due to the enhanced 
9 grid, the resolution in the far-field jet region is ~ 10 
times larger. Also, with the use of a parabolic interpola- 
tion scheme, the overall effective resolution is additionally 
enhanced. Compared to our previous model this gives us an 
effective 6 resolution of « 9000 zones. 

All the results in this paper are robust to changes in 
numerical resolution, as has been explicitly verified by con- 
vergence testing in which the resolution, order of the spa- 
tial interpolation, and #-grid geometry were varied. This in- 
cluded varying the transition radius at which the grid reso- 
lution in the jet becomes larger than the disk via different 
Qj = {1.3, 1.8}, using resolutions of 256 x 256 and 512 x 256, 
and using linear and parabolic spatial interpolation schemes. 
In all cases, no qualitative differences were found. 



4.2 Model Setup 

All the experiments performed evolve an initially weakly 
magnetized torus around a Kerr black hole in axisymmetry. 
A generic model is used so the results mostly can be ap- 
plied to any black hole system. The model of the initial disk 
is iden tical to the model studied in IMcKinnev fc Gammid 
i2004h . and so direct comparisons can be made to our prior 
work. The new models have an extended radial grid and 
have a longer duration. These were required to study the 
large-scale propagation of the Poynting-dominated jet. 

A black hole spin of j = 0.9375 is chosen, but this pro- 
duces sim ilar results to models with .5 < j < 0.99 (see 
also, e.g., IMcKinnev fc Gammid 120041) . This value of the 
black hole spin is close to the equilibrium black hole spin 
for a black hole accreting a disk with a height (H) to radius 
(R) r atio of H/R ~ 0.26 iGammie. Shapiro, fc McKinnevI 

The model is run for At = 1.4 x 10 4 i 9 , which is about 
52 orbital periods at the pressure maximum and about 1150 
orbital periods at the black hole horizon. 

Since the model is axisymmetric , disk turbulence is not 
susta ined after about t ~ 3000t 9 iMcKinnev fc Gammid 
|2004[) . after which the anti-dynamo theorem prevails 
JCowlindll93l . However, while this affects the disk accre- 
tion by forcing it to become more "laminar," this does not 
affect the evolution of the Poynting-dominated jet. That is, 
from the time of turbulent accretion to "laminar" accretion, 
the funnel region is mostly unchanged. Indeed, the far-field 
jet that has already formed is causally disconnected from 
the region where the accretion disk would still be turbulent. 



4.3 Initial Conditions 

The initial conditions consist of an equilibrium 
torus, which is a "donut " of plasma with a black 
hole at the center iFishbone fc Moncried Il976l : 
lAbramowicz. Jaroszinski. fc SikoraT^TsT) ! This torus is 
a solution to the axisymmetric stationary hydrodynamic 
equations and is supported against gravity by centr ifugal 
and p ressure forces. The solutions of iFishbone fc MoncrieJ 
( 1976) corresponding to = const, are used. The initial 
inner edge of the torus is set at r e dge = 6. The equation of 
state is a gamma-law id eal gas with 7 = 4/3, but o ther 7 
lead to similar results iMcKinnev fc Gammid 120041,) ■ The 
models have v/u$ = 4.281, the pressure maximum is located 
at Tmax = 12r 9 , the inner edge at (r,8) = (6r 9 ,7r/2), and 
the outer edge at (r,9) = (42r 9 ,7r/2). The orbital period at 
the pressure maximum is 2ir(a + (r max /r g f'^ 2 )t a ~ 267t 9 , 
as measured by an observer at infinity. The torus chosen 
has a disk height (H) to radius (R) ratio of H/R ~ 0.26 
on average, and is slightly thinner (thicker) at the inner 
(outer) radial edge. 

As done by IMcKinnev fc Gammid i2004Tl . a purely 
poloidal magnetic field is placed into the initial torus. The 
field can be described using a vector potential with a single 
nonzero component, and we choose A$ oc MAX(po/po,m<i:c — 
0.2,0). The field is therefore restricted to regions with 
po/po,max > 0.2. The restriction of the field to be inside 
the torus is to avoid regions that are numerically estimated 
to have large magnetic field per unit particle just outside the 
torus. The field is normalized so that the minimum ratio of 
gas to magnetic pressure is 100. The equilibrium is therefore 
only weakly perturbed by the magnetic field. 

While the torus is stable to axisymmetric hydrody- 
namic instabilities, the torus is un stable to global non- 
axisy mmetric hydrodynamic modes dPapaloizou fc Pringld 
1983). However, when the torus is embedded in a weak mag- 
netic field, the magnetorotational instability (MRI) dom- 
inates those h ydrodynamic modes fealbus fc Hawlevlll99ll : 
lGammiel2 004). Small perturbations are introduced in the in- 
ternal energy density, which induces velocity perturbations 
that seed the instability. 

No initial large-scale net vertical field is necessary 
to drive a Poynting jet, since a large-scale poloidal field 
is self-consistently generated with t he above prescrip- 
tion for the initial f ield g e ometry ij Hirose^t alJ 120041 : 
IMcKinnev fc Gammid 12004 IMcKinnev! l2005al) . In ad- 
dition, other field geometries, such as a net vertical 
field or multipl e field loops within the d isk, lead to 
similar results (IMcKinnev fc Gammid 120041) . Reconnec- 
tion efficiently erases the i nitial geometrical differences 
iMcKinnev fc Gammid f2004l : lMcKinnevlf2005al) . Only con- 
trived fields wit hout any poloidal com ponent lead to 
no Poynting jet jDe Villiers et all l2005ah . Otherwise, the 
magneto-rotational instability drives the amplification of 
toroidal and poloidal flux through differential rotation. 

Beyond the disk, we choose two different models of the 
matter. The first is to set the densities to the "floor model" 
value (denoted as environment 1). The second is to set the 
rest-mass density to po = 10~ 4 r _3//2 po,disfc and the internal 
energy density tou= 10 _6 r _5 ^ 2 po,disk (denoted as environ- 
ment 2). In both cases, there is no magnetic field outside the 
disk. The velocity everywhere is set to be equal to the free- 
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falling frame. The fiducial model described below focuses on 
the second model of the surrounding medium, while the first 
model is used to demonstrate the effect of the existence of 
an environment on the jet. 

4.4 Boundary Conditions 

Two different models are chosen for the outer bound- 
ary condition. The first model is to use a so-called "out- 
flow" boundary condition, for which all primitive variables 
({po,u, u l , B 1 }) are projected into the ghost zones while for- 
bidding inflow (associated with environment f). The second 
model is to inject matter at some specified rate at the outer 
boundary, unless there is outflow. For the second model, 
unless there is outflow the outer boundary is set to inject 
mass at the free-fall rate with no angular momentum, with 
a density of p = f 0~ 4 r~ 3/2 po,disfc and u = f 0~ 6 r~ 5/2 p , diak 
(associated with environment 2). 

The fiducial model below uses the second model, while 
the first model is used together with the "floor model" ver- 
sion of the initial conditions in order to investigate the effect 
of the environment on the jet. 

4.5 Astrophysical Applicability of Model 

As applied to the collapsar model, the outer grid radius of 
r ou t — 10 4 r o corresponds to about 20 presupernova core 
radii ijWooslev fc Weaverlf 99sh . ~ f0 10 km, or about f/fOth 
the entire star's radius. As applied to M87, the outer radius 
of the grid is at 1.4pc. For x-ray binaries, the outer radius 
of the grid is at O.OOf AU. 

For the collapsar model, the chosen black hole spin 
would be achieved in the l ate phase of jet production 
jMacFadven fc Woosle"vlll999l) . As applied to AGN and x- 
ray binaries, the results should be approximately applicable 
to systems with 0.5 < j < 0.99. 

The duration of the evolution for the collapsar model is 
~ 0.2 seconds, and at the spin chosen would correspond to 
the late phase of accr etion when the Poynting flux reaches 
its largest magnitude (fMacFadven fc Wooslevll999l) . For the 
AGN M87 the duration corresponds to ~ 7 years. For the 
x-ray binary GRS f9f5+f05 the duration corresponds to 
~ f second. 

The chosen H/R and the extent of the disk are similar 
to the acc retion disk that is expected to form in the collap- 
sar model jMacF adven fc Wooslevll999llKohri fc Mineshigel 
120021: iKohri et alJl2005h . While in AGN and x-ray binary 
systems an extended disk should be studied, the extent 
of the disk does not affect the results of the jet pro- 
duced since the time scale of the simulations are relatively 
short compared to the outer disks' accretion time scale, 
t ~ 2ty j \olQ,k) ~ 10 5 t 9 , where a ~ O.Of and Qk is the Kep- 
lerian angular frequency. The disk thickness is similar to the 
disk thickness of radiatively inefficient disk models that are 
used to represent the disks in some AGN a nd x-ray bina- 
ries d uring their radiatively inefficient cycle dNaravan fc Yil 
If 9951) . A study of the effect of disk thickness is left for future 
work. 

The "Bondi-flow" model of the matter beyond the disk 
(environment 2) can be compared to the collapsar model. 
This would correspond to the collapsing envelope of a mas- 
sive star. Presupernova models suggest that the infalling 



matter is at about 30% the free-fall rate we have chosen 
iMacFadven fc Wooslevll 999), but this small percent differ- 
ence is unlikely to significantly affect the jet formation. In 
addition, the angular momentum of the envelope is not im- 
portant since the disk is already present in our model and 
the timescale for adding mass to the disk in an accretion 
shock is much longer than the duration of the simulation. In 
collapsar models, the density structure near the equatorial 
plane varies between po oc r~ 3 ^ 2 to pp oc r~ 2 and the inter- 
nal en ergy is u oc r~ 5 ^ 2 to u oc r~ 2 ' 7 jMacFadven fc Wooslevl 
f999). This is similar to our model, which is expected since 
the collapsar model involves mostly spherical collapse. As 
applied to AGN and x-ray binaries, the second model (en- 
vironment 2) is a type of Bondi inflow of the surrounding 
medium. 

The fiducial model below uses the second model, while 
the first model is used together with the "floor model" ver- 
sion of the initial conditions in order to investigate the effect 
of the environment on the jet. 



5 NUMERICAL RESULTS 

The overall character of the accret ion flow is u nchanged 
compa red to the descriptions given in lMcKinnev fc Gammiei 
(2004). The disk enters a long, quasi-steady phase in which 
the accretion rates of rest-mass, angular momentum, and 
energy onto the black hole fluctuate around a well-defined 
mean. Beyond t ~ 3000i 9 , the mass accretion rate drops, 
but all ot her quantities in the disk and jet remain similar. 

As in lMcKinnev fc Gammiei i2004h . a mildly relativistic 
wind is launched from the inner edge of the disk with a 
half-opening angle of roughly f6° to 45°. In this model the 
disk wind has T ~ 1.5. However, a long temporal study of 
this disk wind depends on sustaining disk turbulence, which 
decays beyond t ~ 3000t 9 in this axisymmetric model. 

The coronal-funnel boundary contains shocks with a 
sonic Mach number of M s ~ 100. The inner-radial interface 
between the disk and corona is a site of vigorous reconnec- 
tion due to the magnetic buoyancy and convective instabil- 
ities present there. These two parts of the corona (coronal- 
funnel wall and inner-radial disk-corona interface) are about 
100 times hotter than the bulk of the disk. Thus, these coro- 
nal components are a possible sites for Comptonization and 
nonthermal particle acce leration. 

Meanwhile, as in iMcKinnev fc Gammiei J2004h . a 
Poynting-dominated jet has formed. The Poynting- 
dominated jet forms as the differential rotation of the disk 
and the frame-dragging of the black hole induce a signif- 
icant toroidal field that launches material away from the 
black hole. The Poynting-dominated jet is launched once 
the pressure of the mostly nonmagnetic initial funnel mate- 
rial is lower than the toroidal magnetic pressure. This occurs 
within t < 500t g . 

Figure and figure [5] show the log of rest-mass den- 
sity and magnetic field at the final time of t = 1.4 x 10 4 i 9 . 
Clearly, the jet has pummelled its way through the surround- 
ing medium. By the end of the simulation, the field has 
been self-consistently launched into the funnel region and 
has a regular geometry there. The density and field show 
how the surrounding medium is shocked by the jet in a lat- 
eral direction. In the disk and at the surface of the disk the 
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Figure 1. Panel (A) shows final distribution of logpo on the 
Cartesian plane. Black hole is located at center. Panel (B) 
shows magnetic field overlaid on top of log of density. Outer 
scale is r = 10 4 r 9 . Time is t = 1.4 X 10 4 t 9 . Color represents 
log(po/po disk) with dark red highest and dark blue lowest. The 
final state has a density maximum of po ps 2po ^uk an d a min- 
imum of po ~ 10 -13 po !( jisfc at large radii. Grid zones are not 
smoothed to show grid structure. Outer-radial zones are large, 
but outer 6 zones are below the resolution of the figure. For the 
purpo ses of properly visualizi n g the accretion flow and jet, we 
follow MacFadvcn Sz Wooslcv 1 1999) and show both the nega- 
tive and positive x-region by duplicating the axisymmetric result 
across the vertical axis. 
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Figure 2. Same as figure HI but outer scale is r = 10 2 r 9 . 



field is curved on the scale of the disk scale height. Within 
r < 10 2 r 9 the funnel field is ordered and stable due to the 
poloidal field dominance. However, beyond r ~ 10 — 10 2 r 9 
the poloidal field is relatively weak compared to the toroidal 
field and the field lines bend and oscillate erratically due 
to pinch instabilities. The radial scale of the oscillations is 
10 2 r 9 (but up to 10 3 r 9 and as small as 10r 9 ), where r ~ 10r 9 
is the radius where poloidal and toroidal field strengths are 
equal. By the end of the simulation, the jet has only fully 
evolved to a state independent of the initial conditions out 
to r ~ 5 x 10 3 r 9 , beyond which the jet features are a result 
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Figure 3. Contours for the disk-corona, corona-wind, and wind- 
jet boundaries at t fa 1500t g and an outer scale of r rj 10 2 r 9 . 
The disk-corona boundary is a cyan contour where /3 ~ u/b 2 = 3, 
the corona- wind boundary is a magenta contour where (3 = 1, 
and the wind-jet boundary is a red contour where b 2 /(poc 2 ) = 1. 
The black contour denotes the boundary beyond which material 
is unbound and outbound (wind + jet). 



of the tail-end of the initial launch of the field. The head of 
the jet has passed beyond the outer boundary of r = 10 4 r 9 . 

Figure [5] shows that the strongest magnetic field is near 
the black hole in an X-configuration due to the Blandford- 
Znajek effect having a power dPj et /dO oc sin 2 8 that is trun- 
cated by the disk near the equatorial plane. 

Figures [3] and [I] show the energy structure of the disk, 
corona, an d jet. This figure is comparabl e to the left panel of 
figure 2 in iMcKinnev fc Gammia i2004f) . These plots show 
the boundaries between the "disk," "corona," "disk wind," 
and "Poynting jet." The disk-corona boundary is defined 
to be where f3 = p g /pb = 2(7 — l)u/b 2 = 3, the corona- 
wind boundary is defined to be where /3 = 1, and the wind- 
jet boundary is defined to be where b 2 /(poc 2 ) = 1. The 
plot also shows the boundary beyond which material is un- 
bound and outbound (wind and jet), which here is defined 
as where u p > and u t = — 1. This corresponds to where 
non-interacting particles are unbound, although thermal and 
magnetic energy can further unbind the matter. Hence, more 
material within the surface is also unbound. 

The jet stays collimated, while the disk wind has a large 
opening angle at large radii. Beyond r ss 10r 9 , the jet under- 
goes poloidal oscillations due to toroidal pinch instabilities. 
By r > 100r s , pinch instabilities subside once magnetic en- 
ergy is converted into thermal energy, which then supports 
the jet. Residual slow dense and fast diffuse patches from 
the instability are present in the jet, such as the slower and 
cooler dense blob shown at top left corner of figure. 

Figure |S] shows the disk, corona, and jet magnetic field 
structure during the turbulent phase of accretion at t « 
1500£ s . Compared to figure this shows the turbulence in 
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Figure 4. Same as figure l3l but for t Ri 1.4 X 10 4 c 9 and an 
outer scale of r 53 10 3 r 9 . At large scales, the cyan and magenta 
contours closer to the equatorial plane are not expected to cleanly 
distinguish any particular structure. 
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Figure 5. Field geometry near black hole for t as 1500i fl during 
phase of strong disk turbulence. 



the disk, but is otherwise similar. The jet, disk, and coronal 
structures remain mostly unchanged at late times despite 
the decay of disk turbulence. That is, the current sheets 
in the disk do not decay and continue to support the field 
around the black hole that leads to the Blandford-Znajek 
effect. The coronal thickness and radial extent do not require 
the disk turbulence, which only adds to the time-dependence 
of the disk wind. Note that despite the fact that the disk 
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Figure 6. Plot of log Too with red highest (Too ~ 10 4 ) and blue 
lowest. Yellow is Too ~ 10 2 — 10 3 . Panel (A) shows outer scale 
of r = 10 4 r 9 . Panel (B) shows outer scale of r = 10 3 r 9 with 
same color scale. Jet is independent of initial conditions by r fa 
5 X 10 3 r 9 . Jet becomes conical at large radii with a core half- 
opening angle 9j pa 5° . 



wind only contains small-scale magn etic fields, such fields 
still lead to self-collimation ( 020021) . 

Figure|S|shows Too, which reaches up to Too ~ 10 3 -10 4 , 
for an outer scale of r = 10 4 r 9 (panel A) and an outer 
scale of r = 10 3 r 9 (panel B). The inner-radial region is not 
shown since Too is divergent near the injection region where 
the ideal-MHD approximation breaks down. Different real- 
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Figure 7. A radial cross section along a mid-level field line in the 
jet showing the velocity structure. The top panel shows T^ M ' 
(solid line), T^f A ' (dotted line), and T (dashed line). The mid- 
dle panel shows u^ 1 ^ 1 (solid line), (dotted line), and 

(dashed line). The bottom panel shows |u r | (solid line) and u* > 
(dotted line). 

izations (random seed of perturbations in disk) lead up to 
about Too ~ 10 4 as shown for the lower pole in the color fig- 
ures. This particular model was chosen for presentation for 
its diversity between the two polar axes. The upper polar 
axis is fairly well structured, while the lower polar axis has 
undergone an atypically strong magnetic pinch instability. 
Various realizations show that the upper polar axis behav- 
ior is typical, so this is studied in detail below. The strong 
hollow-cone structure of the lower jet is due to the strongest 
field being located at the interface between the jet and the 
surrounding medium, and this is related to the fact that the 
BZ-flux is oc sin 2 6, which vanishes identically along the po- 
lar axis. It is only the disk+corona that has truncated the 
energy extracted, otherwise the peak power would be at the 
equator. As described in the next section, at larger radii the 
jet undergoes collimation that results in more of the energy 
near the polar axis. 

5.1 Radial Jet Structure 

Figure shows the velocity structure of the Poynting- 
dominated jet along a mid-level field line, which starts at 
6j « 46° on the black hole horizon and goes to large dis- 
tance, where the funnel region starts at 6j « 63° on the 
horizon. The quantities shown are time-averaged values at 
late-time. 

The top panel of figure Q shows the electromagnetic 
(EM) energy flux per unit rest-mass flux (r£^ ), the mat- 
ter (MA) energy flux per unit rest-mass flux (r^ IA '), and 
the Lorentz factor as measured by an observer at infinity 
(r). For 2 < r < 10, the value of F^ M ' is highly oscillatory. 
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This shows the location of the stagnation surface, where the 
poloidal velocity u p = 0, is temporally variable. The stag- 
nation surface is, at least, where the ideal-MHD approxima- 
tion breaks down and thus at any moment the value of F^? M ' 
nearly diverges. Notice that while at r < 10 3 r 9 the Poynting 
flux dominates, the Poynting flux energy is slowly converted 
into enthalpy flux due to nonlinear time-dependent shock 
heating. Shock s are expected beyond the ma gneto-fast sur- 
face (see, e.g., iBoeovalov fc Tsinganosl l2005f) and here are 
partially due to pinch in stabilities feichleilll993l : iBegelmanl 
ll998l:ISikora et alJl2005f) . 

For r > 10 3 r 9 , the enthalpy flux and Poynting flux are 



in equipartition (r 



(MA) 



i(BM) 



) and so an equipartition 



"magnetic fireball" has formed. The terminal Lorentz factor 
is of order Too ~ 400 for this choice of flow line. The Lorentz 
factor by r ~ 10 4 is still only 3 < T < f 5, with smaller spa- 
tial patches with T ~ 25. For a radiation pressure-dominated 
and optically thick flow, this implies there will be an ex- 
tended acceleration region where the magnetic fireball loses 
energy during adiabatic expansion. 

In a narrow region within the core of the jet (not shown 
in figure within a half-opening angle of 6 ~ 5°, the jet 
accelerates with 

F~( I ^)°- 3 , (inner) (22) 

for 10 < r < 5 x 10 3 r 9 , after which the initial conditions 
dominate the flow. Thus r ~ 6 has been reached in the 
core. The middle and outer angular parts of the jet follow 

0.45 

r ~ ( — I , (inner) (23) 



for 3 < r < 10 3 r 9 , after which the field becomes pseudo- 
conical and there is no sustained acceleration (or the flow 
even decelerates). Thus, F ~ 10 has been reached. As dis- 
cussed below, the field is nearly paraboloidal for the radial 
range over which there is such power-law acceleration. The 
behavior of T(r) is in good agreement with an analytic study 
of acceleration in a parab oloidal field that predicts F tx r 1 ^ 2 
jBeskin fc Nokhrinal2005l) . Thus our numerical result is con- 
sistent with their analytical result. 

The next lower panel of figure [7| shows that the electro- 
magnetic angular momentum is converted to enthalpy an- 
gular momentum and that there is a conversion to particle 
angular momentum. 

The bottom panel of figure |7| shows the approximate 
orthonormal radial (u r « yjg rr u r ) and <f> (ir « ^/g^^u^) 
4- velocities in Boyer-Lindquist coordinates, where the ap- 
proximation being made is that the metric is dominated 
by the diagonal terms for r > 10r 9 . The radially directed 
motion becomes relativistic, while the (^-component of the 
4-velocity becomes sub-relativistic at large distances. How- 
ever, the fluid and field rotation is still large enough that 
it may stabilize the jet against m = 1 kink instabilities, 
as shown in the nonrel ativistic case l|Ouved et al.l 2003; 
iNakamu ra fc MeieJl200 4[) an d analytically in the relativistic 
case llTomhri£feu^t^ni200 if) . 

Figure |H| shows the collimation angle, rest-mass and in- 
ternal energy densities, ratio of magnetic to rest-mass en- 
ergy density, toroidal field strength in orthonormal Boyer- 
Lindquist coordinates, and the pitch angle along a mid- level 
field line in the Poynting-dominated jet. The field line is 



the same as that in figure |7| and now for both poles of the 
simulation. 

The two poles are similar but not necessarily identi- 
cal. Many simulations that were performed show at least 
one polar jet develops some pinch instabilities. In the model 
shown, one side develops very strong instabilities by the end 
of the simulation. The development of the pinch instability 
mostly affects the far-radial Lorentz factor (F and Too - 
and so densities). The value of V is up to 5 times larger in 
pinched regions than non-pinched regions, while Too can be 
up to a factor 10 times larger in pinched regions compared 
to non-pinched regions. 

In figure 

HI for 7r 9 < r < 100r 9 there is a region of 
collimation slightly faster than power-law. For the field line 
closest to the disk wind starting at 9j w 57°, collimation is 
nearly a power-law with 



57° 



2.8r c 



(24) 



8° beyond. Closer to the polar 
axis the collimation is faster. For the field line starting at 
Oi « 17°, 

6j oc r-° A , (25) 

approximately up to r ~ 120r 9 and 9j ~ 5° beyond. The in- 
ner core for 5 < r < 10 3 r 9 follows 0j oc r -1 / 2 , a paraboloidal 
field. See also discussion related to figure [TTjl 

From a measurement of the forces along and across the 
field lines within the Poynting jet, the inner-radial collima- 
tion is dominated by confinement by the disk+corona, while 
in the mid-radial range the disk wind collimates the Poynt- 
ing flow. Far from the disk wind, the internal jet is collimated 
by internal hoop stresses due to the paraboloidal field. Note 
that the classical hoop-stress paradigm that jets can self- 
collimate is not fully tested here, but these results suggest 
that collimation is in large part due to, or at least requires, 
the disk and disk wind. 

Notice from figure |7| that there is little acceleration be- 
yond r ~ 10 3 r 9 . It is well known t hat magnetic ac celera- 
tion requires field lines to collimate llEichlerl Il99~3l : Oil 9931 : 
iBegelman &Tilll994l: | Tomimatsulll994l) . In a conical flow, 
the transfer of magnetic energy into kinetic energy is in- 
efficient. As shown in figure |U after r > 10 3 r 9 the flow 
oscillates around a pseudo-conical asymptote with a mean 
half-opening angle of Oj ~ 5°. Beyond this radius, the mag- 
netic energy is instead converted to thermal energy. The 
pseudo-conical asymptote results once the magnetic energy 
no longer completely dominates the other energy scales and 
the disk wind no longer helps collimate the flow. The final 
opening angle should be similar for other black hole spins, 
as found in prior studies jMcKinnev fc Gammieil2004D . 

The lower panels of figure |H| show the toroidal field and 
pitch angle (a = tan _1 (|B r |/|i?'' > |)), which shows that even- 
tually the toroidal field completely dominates the poloidal 
field, and the toroidal field remains ordered. The toroidal 
field within r < 390r 9 is well fit by 



diskC 



= [G] 



-0.0023 



3907-, 



(26) 



and is otherwise well fit by 
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Figure 8. A radial cross-section of both poles (solid and dotted lines) of the jet along a mid-level field line in the jet showing the jet 
opening angle in degrees, density, and magnetic structure (B r and \B^\ in Gaussian units). In the density plots and the plot of B^, the 
overlapping dashed lines are power-law fits. The pitch angle (lower right panel) in Boyer-Lindquist coordinates (which vanishes at the 
horizon) has an overlapping fit (dash-dotted line) to the Tomi matsu et all l200ll) kink stability criteria when including field rotation. 



B" 



diskC 



-AG] 



-0.0023 



390r o 



(27) 



The transition radius is 



390r o . This transition is as- 



sociated with the field line going from collimating to not 
collimating (or slightly decollimating) . This is along one 
mid-level field line, and the coefficients and power laws are 
slightly different for each field line. 

The pitch angle shows that the field becomes toroidally 
dominated and at large radius the magnetic loops have a 
pitch angle of < 1° by 10 3 r 3 . The small pitch-angle sug- 
gests that the flow may be kink unstable. The Kruskal- 



Shafranov criterion (|S <#> /B I '| > sj 94,4, /gee), where ^Jgee is 
approximately the vertical length of the jet, is the standard 
kink instability criterion. However, a self-consistent treat- 
ment determined that a Poynting-dominated jet is kink un- 
stable only if both the Kruskal-Shafranov criterion is satis- 
fied and for a mostly rad ial flow that \B' t '/B r \ > ^/glp^/Rh 
llTomimatsu et aT]l200ll) . where Rl = c/Qf is the light 
cylinder radius and cylindrical radius R ~ ^fg~H> at large 
distances. Their result only strictly applies inside the light 
cylinder and for a uniform field, but their result is sugges- 
tive of the general result. The lower right panel of figure |H1 
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shows tan - (Rl/ y/g<t,</>) overlapping the pitch angle, and so 
the flow is marginally stable to the kink instability. Thus, 
the jet is not expected to be (violently) kink unstable if sim- 
ulated in 3D. 

The above discussion implies that the radial magnetic 
field is well fit by 



B r 



-B 4 



(28) 



where we also find that flp w fi# /2 is approximately con- 
stant along a field line (also see, e.g.. lMcKinnev fc Gamniiel 
|2004) . where Qh = j/(2r_fi) is the angular frequency of 
the black hole. The marginal kink stability of our solu- 
tion is also in concorda nce with the fact that fif ~ fin/2 
llTomimatsu et all200lfl . since a Poynting-dominated jet fol- 
lowing such a rotation was shown to be marginally kink sta- 
ble. 

As pinch instability-driven oscillations develop, this 
leads to arbitrarily sized patches moving at arbitrary rel- 
ative velocities as required by any internal shock model. A 
patch is simply defined as a localized increase or decrease 
in various fluid quantities. For example, figure |H] shows that 
at large radii that the percent difference in the fluctuating 
part of the opening angle depends on the hemisphere and 
varies between 10% to 50%. Comparing against each quan- 
tity smoothed over a logarithmic interval of dlog(r) = 0.4 
for the hemisphere with less oscillations, the typical percent 
difference of quantities along the length of the jet along a 
mid-level field line is: 10% for po, b 2 , B r , B^, T, and u$ ; 
20 — 50% for u ; and 30% for Too. The other hemisphere with 
larger 8j oscillations has w 5 times the percent difference 
compared to the first hemisphere. For both hemispheres, 
the largest several patches have an order unity percent dif- 
ference. The typical size of a patch is ~ 100r 9 to ~ 10 3 r g in 
the lab frame along the length of the jet. 

The top two right panels of figure |H] show the rest-mass 
density per disk density and internal energy density per rest- 
mass density. The middle left panel shows the ratio of (twice) 
the comoving magnetic energy density per unit rest-mass 
density. The region within r < 6r 3 is dominated by the 
floor model, and this is where the mass is injected into the 
Poynting-jet. This explains the kink in the density profiles 
at r ~ 6r 9 and the plateau in b 2 /po- For r > Wr g , the floor 
model is never activated and the matter in the Poynting 
jet is self-consistently evolved. The inner-radial rest-mass 
density is well fit by 



P0 ~ 1.5 x 10- 9 



PO.disk 



120r c 



(inner) 



(29) 



The outer-radial rest-mass density is well fit by 

-2.2 

P0 ~ 1.5 x 10" 9 [ -rZr- I • (outer) (30) 



PO.disk 



120r £ 



internal energy density is moderately fit by 

-1.8 

(inner) (31) 



JL^ « 4 .5 x 10" 9 ' T 



120r c 



The outer-radial internal energy density is moderately fit by 
4.5 x 10" 9 ( -rZr- I • (outer) (32) 



PO.diskC' 



120r c 



For both the rest-mass and internal energy densities, 
the transition radius is r ~ 120r s , which is near the tran- 
sition to a superfast flow. The location of this transition 
may be due to the thickness of the disk, since the disk leads 
to collimation at small radii. A flow that co llimates more 
leads to the fast surface closer to the o rigin feichledll993l : 
iBeeelman fc Lilll994t lTomimatsulll994h . Hence, if a thick 
disk is required for collimation, then a thick disk is required 
for a superfast transition at small radii. 

The internal energy density stabilizes in value once the 
equipartition "magnetic fireball" has formed by r ~ 10 3 r 9 . 
However, the simulation did not extend far enough to guar- 
antee that the equipartition extends for much larger radii. 
The internal energy may continue to grow, stay in equiparti- 
tion, or decrease at the expense of adiabatic expansion and 
acceleration of the jet. The resolution if this issue is left for 
future work. 



5.2 9 Jet Structure 

Figure[5]shows the 9 cross-section of the jet at r = 5 x 10 3 r s 
for the upper polar axis according to figuresQ |5] and|5] This 
is the location by which the jet has stabilized in time, where 
farther regions are still dependent on the initial conditions. 
The hot fast "core" of the jet includes 9 < 5° at this radius 
and is marked by the left dashed vertical line. Also, note 
that 9 « 8° corresponds to the radius obtained for the "mid- 
level" field line shown for radial cross-sections in figures |7||S] 
and [Hll 

The right dashed vertical line marks the boundary be- 
tween the last field line that connects to the black hole. This 
extended angular distribution out to 9 < 27° is colder than 
the hot core of the jet. Beyond this region is the surrounding 
infalling medium. The disk wind has not reached such large 
radii. 

The very inner core of the jet is denser and slower. This 
feature is simply a result of the fact that the Blandford- 
Znajek flux is tx sin 2 9. Most of the Poynting energy flux is 
at the interface between the jet and surrounding medium at 
small radii, where at larger radii the jet internally evolves to 
collimate into a narrower inner core where the energy flux 
is concentrated. 

In most of our simulated models, the jet structure is 
nearly Gaussian. The fiducial model has 



e(9) 



e e 



-e 2 - /2ej 



(33) 



where eo ~ 0.18Mc 2 and 9 t ~ 8°. The total luminosity per 
pole is Lj « 0.023M c 2 , where 10% of that is in the "core" 
within a half-opening angle of 5° with the largest Lorentz 
factor within the jet. Also, Too is approximately Gaussian 
with 



1^(0) = roo,oe 



-e 2 /2ej 



(34) 



where Too,!) ~ 3 x 10 and 9^ « 4.3°. The Lorentz factor is 
approximately Gaussian with 

r(6») = r oe - 02/2B r, (35) 

where To « 5 and 9r ~ 11°. 

The jet structure can also be modelled by a small uni- 
form core with sharp exponential wings. A good fit is 



e(9) = O.lMc 2 (9 s: 11°) 



(36) 
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Figure 9. A 9 cross-section of the jet at r = 5 X 10 3 r 9 showing the velocity, density, energy, and magnetic structure in Gauss. Too and 
t(9)/Mc 2 have overlapping Gaussian fits (dashed lines) and overlapping uniform+exponential fits (dotted lines), where e(9) = —r 2 T[ 
r 2 pof r TcxD (Boyer-Lindquist coordinates). A Gaussian fit to T is shown as an overlapping dotted line. Lower right corner shows e/Mc 2 
with overlapping power-law fit (dashed line). The magnetic field components are in Boyer-Lindquist coordinates in an orthonormal basis. 
Vertical dashed lines marks the outer jet and jet core. Those panels with two separate y-axes labels have the labels vertically ordered by 
the function value near 9 = 0. 



and otherwise 



e(0) = 0.1M 



c 2 e (-3(9/S e -l)) 



(<9 ^ ll c 



where 6, = 11°. We find that 
T x = 3 x 10 3 (6 <: 5.8°) 
and otherwise 

r^3x ioV-W-- 1 " (0 5.8°), 
where Ooo = 5.8°. 



(37) 
(38) 
(39) 



The bottom right corner more clearly demonstrates the 
fact that there is a significant amount of energy flux at low 
Lorentz factors. The dependence roughly follows 



Mc 2 



Toe 
103 



1/5 



(40) 
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5.3 Jet Characteristic Structure 

In this section the characteristic structure of the jet is con- 
sidered . First, the ide al MHD dispersion relation is given, 
e.g., in lGammie et alJ <l2003al) (there is a sign typo there), 
and summarized here. In the comoving frame, the dispersion 
relation is 

= £>0 M ) 

= lu (lu 2 - k 2 a ) 

x (lu 4 - lu 2 (K 2 c 2 ms + c s k 2 a /c) + K 2 c 2 s kl) , (41) 

where k a = (k • v A ), c 2 ms = (\\ + c 2 (l - v A /c 2 )) is the 
magnetosonic speed, c 2 = (d(p + uj/dp)^ 1 is the relativistic 
sound speed, va = B/\/£ is the relativistic Alfven veloc- 
ity, £ = b 2 + w, and w = p + u + p. The invariant scalars 
defining the comoving dispersion relation are lu — —k^u^, 
K 2 = Ky.K" = fc M fc" + lu 2 , where K„ = P^k" = k„ - luu^ 
is the projected wave vector normal to the fluid 4-velocity, 
v| = b M & M /£, and (k • v A ) = k^/VE. The terms in the 
dispersion relation correspond to, respectively from left to 
right, the zero frequency entropy mode, the left and right go- 
ing Alfven modes, and the left and right going fast and slow 
modes. The eighth mode is eliminated by the no-monopoles 
constraint. The dispersion relation gives the ingoing and out- 
going slow, Alfven and fast surfaces. 

Energy can be extracted from the black hole if 
and only if the Alfv en point lies inside the ergosphere 
(Takahashi ct al Optimal acceleration of the flow by 

conversion of Poynting flux to kinetic energy flux occurs be- 
yond the outer fast su rface iEichlerl Il993t iBeeelman fc Lil 
11994 lTomimatsulll994h . Other surfaces of interest include: 
the horizon at rg = 1 + \J\ — j 2 ; the ergosphere at 
r = 1 + yl — (j cos 9) 2 ; the coordinate basis light surface in 
where ^fg&t, = c/Qf, where asymptotically ^/g^p = rsin(0) 
is the Minkowski cylindrical radius; and the surface in Boyer- 
Lindquist coordinates where the toroidal field equals the 
poloidal field (\B^\ = \B r \). Finally, there is a stagnation 
surface where the poloidal velocity u p = 0. In a field con- 
fined jet where no matter can cross field lines into the jet, 
and if the jet has inflow near the black hole and outflow far 
from the black hole, then this necessarily marks at least one 
location where rest-mass must be created either by charge 
starving the magnetosph ere till the Goldreich - Julian charge 
density is reached jGoldreich fc Julian! lilii) . or pair pro- 
duction rates sustains the rest-mass density to a larger value 
than the Goldreich- Julian value. 

Figure I1UI shows 6j as a function of radius for a field 
line close to the polar axis, in the middle of the jet, and 
one close to the outer edge of the jet. Also shown is the 
Blandford-Zna jek parabolo i dal so lution given by equation 
7.1 in iBlandford fc Znaiekl lll977li . The three paraboloidal 
field lines correspond to X(r,8)/C = 1.6,1.2,0.75 in that 
equation, where X(r,9)/C = 2(1 — log(2)) along the po- 
lar axis. Note that at large radii that a paraboloidal field 
has 6j oc r -1 / 2 , cylindrical has 6j tx r _1 , and conical has 
9j oc r°. The ingoing fast point is nearly coincident with the 
horizon. The stagnation points have time-dependent posi- 
tions that vary within approximately ±30% differences. 

Notice from figure I1UI that the field near the hori- 
zon is nearly monopolar. After the stagnation point the 
field line becomes nearly paraboloidal up to the outgoing 




r c a /GM 

Figure 10. A radial cross-section of the Poynting-dominated 
jet along three field lines showing the opening angle in degrees. 
The field lines are close to the polar axis, along a mid-level field 
line, and close to the outer edge of the Poynting-dominated jet. 
Overlapping dashed lines are paraboloidal fits within the ra- 
dial range that has a similar collimation. Solid squares mark 
the ingoing/outgoing fast points. Solid triangles mark the ingo- 
ing/outgoing Alfven points. Pluses mark the stagnation point. 
Open squares mark where poloidal and toroidal fields are equal. 
Open circles mark the light surface. Dotted lines mark the hori- 
zon and ergosphere. The field lines at larger 8 are shown for a 
shorter radius due to limitations of a program used to extract a 
field line contour. 



Alfven point. Beyond this point, the field is unstable to pinch 
modes and some of the magnetic energy is converted to ther- 
mal energy. For the outer and mid- level field lines, after the 
flow passes through the outgoing fast point, the flow goes 
from nearly paraboloidal to nearly conical. The inner-level 
field line continues to follow a fairly paraboloidal collimation 
with mild oscillations. The magnetic surfaces are therefore 
paraboloidal in the core of the jet and surrounded by conical 
surfaces, which is quite simil ar to the assumed st r uctur e of 
the analytical model bv lTomimatsu &~ Takahashi (l2003ll . 

Figure 1111 shows the characteristic structure of the jet 
for one polar axis. The Alfven surface li es inside the er- 
gosph ere, as required to extract energy jTakahashi et al] 
1990). Clearly the field lines follow nearly a power-law until 
r ~ 10 2 r 9 . The stagnation surface is time-dependent but sta- 
ble. Clearly the transition to a supercritical (superfast) flow 
has occurred. After the fast surface, the field lines stretch 
out and oscillate around a conical asymptote. The fast sur- 
face near the polar axis is at r ~ 250r 9 , where for the other 
hemisphere it is at r ~ 500r 9 . Other features are quantita- 
tively similar. 

Despite the unsteady nature of the flow, the charac- 
teristic structure is quite simple and relatively smooth in 
appearance. This indicates that the flow is mostly station- 
ary. The region near r ~ 10 4 r 9 is still determined by the 
initial conditions, which explains the distorted appearance 
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GRMHD case, we find that the confinement is mostly 
magnetic 1 . 

The radial structure of the jet is also unaffected by the 
environment. For example, one can compare the radial struc- 
ture of the density in the jet (equations 1291 and 1301 to the 
density of the evacuated environment model of the surround- 
ings, 



P0,disk 



= 2 x 10" 



120r c 



(Environment!) (42) 



or to the pre-existing spherical infall environment model of 
the surroundings, 



PO.disk 



8 x 10" 



120r £ 



(Environment2) (43) 



Even with the drastic difference in rest-mass and internal 
energy densities, all the results of the Poynting-dominated 
jet are quantitatively similar. This is much different than the 
results from hydrodynamic simulations mentioned above. 



Figure 11. Poynting-dominated jet characteristic (and other) 
surfaces. Shows a log-log plot of one hemisphere of the time- 
averaged flow at late time. In such a log-log plot, 45° lines cor- 
respond to lines of constant 8. Lines of constant spherical po- 
lar r are horizontal near the 2-axis and vertical near the R- 
axis. The field lines are shown as red lines. From r = out- 
wards: Blue#l: horizon + ingoing-fast ; Cyan#l: ingoing-Alfvcn ; 
Black#l: = \B f \ ; Green#l: ergosphere ; Purple#l: ingoing- 
slow ; Black#:2: stagnation surface where poloidal velocity u p = 
; Purple#2: outgoing-slow ; Cyan#2: outgoing-Alfven ; Black#3: 
\B^\ = \B r \ again ; Green#2: light cylinder ; Blue#2: outgoing- 
fast. The disk and coronal regions have been truncated with a 
power-law cutoff for r < 100r 9 and a conical cutoff for larger 
radii. Within r < 10r 9 the plotting cutoff creates the appearance 
that the fast surface and other lines terminate along the cutoff. 

of the field lines and other artifacts related to the flow being 
unsteady. 



5.4 Dependence on Surroundings 

The jet structure that emerges from the simulation is pri- 
marily due to the internal evolution of the jet rather than 
by the interaction with the surrounding medium. In par- 
ticular, these GRMHD-based results are insensitive to the 
two different models of the initial surrounding medium. One 
model is a surrounding infall of material and the other is 
an "evacuated" exterior region, as discussed in sections 14.31 
and 14.41 The Poynting-dominated jet structure is negligibly 
broader or narrower in the evacuated case due to magnetic 
confinement. The disk wind itself also easily plows through 
the exterior region. 

This is in stark contrast to the simu l ations 
of relativistic hydrodynamic jets (|AJov^j^iL|_|200(i 
[ Zhang. Wooslev. fc MacFadvenl 120031: IZhang W. et alJ 
2004), where they find that the environment is primarily 
responsible for collimating the hydrodynamic flow and 
hydrodynamic collimation shocks help diminish the mixing 
between the jet and the surrounding medium. In the 



6 DISCUSSION 

We find that the Poynting-dominated jet accelerates to 
r ~ 1 along nearly paraboloidal field lines. The maximum 
terminal Lorentz factor is Vac < 10 3 , sufficiently relativistic 
to account for any relativistic jet from GRBs, AGN, or x-ray 
binaries. At large radii, about half of the energy in the jet is 
thermal. Thus, for optically thick flows that are radiation- 
dominated, thermal adiabatic expansion is expected to lead 
to an extended acceleration range far beyond r ~ 10 4 r 9 . 
For optically thin flows that are radiation-dominated, much 
of this energy can be directly released as radiation. In the 
limit that the radiative time-scale is much shorter than the 
jet propagation time-scale, the terminal Lorentz factor is 
limited to ~ V ~ 10. 

The part of the jet with the largest energy flux colli- 
mates from a half-opening angle of 8j ~ 60° near the black 
hole to 9j ~ 5° at large distances. The core of the jet fol- 
lows nearly paraboloidal field lines for all radii and accel- 
erates with approximately F oc r 0,3 from r ~ 10r 9 out to 
r ~ 5 x 10 3 r 9 , after which the initial conditions dominate 
the flow. Thus the core moves with r ~ 6 and may continue 
to accelerate since the field is still collimating by the outer 
radius. The middle and outer angular parts of the jet follow 
nearly paraboloidal lines inside the fast surface and conical 
lines outside the fast surface. The middle and outer parts 
follow r oc r ' 45 from r ~ 3r 9 up to the point where the flow 
becomes conical at r ~ 10 2 r 9 — 10 3 r 9 . Thus these parts move 
with r ~ 10. The acceleration nearly follows the result from 
analytic models wi th exactly paraboloidal fi eld lines that 
found r(r) oc r 1/2 dBeskin fc NokhrinallioO^) . That the ac- 
celeration found is slightly weaker than for paraboloidal lines 
is consistent with the fact that the simulated field lines end 
up not as strongly collimated as paraboloidal. 

1 There is a difference between magnetic confinement and mag- 
netic collimation. Confinement refers to the magnetic field con- 
taining the mass of the jet since there is no hydrodynamic mixing, 
while collimation refers to modifications of the field line geometry 
along the jet. 
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The jet structure is Gaussian-like, so there is a non- 
negligible amount of energy at large angles out to 8j ~ 27°. 
This may explain the observations of some GRBs that ap- 
pear to have a narrow ultrarelativistic component that pro- 
duces the 7-rays and early afterglow and a wide, mildly 
relativistic co mponent that prod uces the radio and opti- 
cal afterglow dBerger et al J 12003). The provided fits to the 
jet structures can be folded into various models, such as 
the probability of obse rving polarised emission in Compton 
drag emission models jGhisellini et al.ll2000l: | Lazzati et alJ 
2004) or whether Compton-drag models work. Consider 
the region within 6 < 27°. This is an expanded, rela- 
tively cold slow portion of the jet. It is possible that the 
gamma-ray burst photons are due to Compton drag from 
soft photons emitted by this jet "sheath" dumping into the 
faster spine feegelman fc Sikorall987tlGhisellini et al.l2000t 
lLazzati et alJl2004h . 



6.1 Toroidal Field Instabilities 

Internal shocks and magnetic instabilities have both been in- 
voked to explain emission from jets associated with GRBs, 
AGN, and x-ray binary systems. We find that the Poynting- 
dominated jet becomes unstable beyond the Alfven surface 
and magnetic energy is converted to thermal energy. Un- 
like during reconnection, this process involves toroidal field 
instabilities that lead to a conversion of magnetic to lat- 
eral kinetic energy. These oscillations drive waves into the 
jet that steepen into shocks and lead to r2f A ' ~ 10 3 . This 
process is consistent with expectations of current-driven in - 
stabilities (lEichlerll993HBegelmanll998l : ISikora et all2005l) . 
These instabilities lead to slightly faster and slower moving 
patches of jet material. Thus, the variability within the jet 
is dominated by toroidal field instabilities. 

Toroidal field instabilities have b een in voked to signifi- 
cantly lower the BZ power output although their 
estimate was based upon the Kruskal-Shafranov criteria for 
kink instability, a nonrelativistic model, and limited the 
BZ power output in the context of an acausal astrophys- 
ical load at large distances. At best their estimate deter- 
mines a limitation of the Poynting flux domination of the 
jet, but not the total jet power or power emitted by the 
black hole. Kink instabilities have also been invoked to de- 
stroy the co herent structure of jets (for a discussion, see, e.g., 
lAppj|l99r3) . and such instabilities are often regarded as suf- 
ficiently strong so that m agnetic confinement is not possible 
( Mill er- Jones et al~l l2006) . However, the Kruskal-Shafranov 
criteria for kink instability is a necessary but not sufficient 
condition for kink instabilities to occur, where rotation of 
the fluid or field lines can stabilize the jet IITomimatsu et alJ 
l200ll) . 

Blandford-Znajek t ype jet solutions natura lly sit near 
marginal kink stability jTomimatsu et aT1 l20011. Since we 
find that the self-consistently simulated Poynting jet only 
marginally satisfies the necessary and sufficient conditions 
for kink instability, the jet is not expected to be (violently) 
kink unstable if simulated in full 3D. This explains why as- 
trophysical jets can go to large distances even if satisfying 
the Kruskal-Shafranov kink criteria. 

Internal shocks are expected to generate some of the 
emission in GRBs, AGN, and x-ray binaries. In the simu- 
lated jet, the patches generated by the toroidal field insta- 



bilities move at large relative Lorentz factors, and so these 
patches can lead to internal shocks. However, by the ra- 
dius simulated, the Lorentz factor is not yet large enough 
for the standard GRB in ternal shock mo d el to be efficient 
jKobavashi et all Il997l : iMeszarosI 120021 iGhirlanda et alJ 
2003; Piran 2005). However, the energy provided by the BZ- 
driven jet may be more than sufficient to allow for inefficient 
shocks, since in section I5~2l we found that L « 0.023A/qc 2 ~ 



4 x 



10 51 erg s" 



for the collapsar model rather than the 
canonical ~ 10 50 erg s _1 for cosmological GRBs. The accel- 
eration region and internal shock region was not simulated 
since the dynamical range required is another ~ 6 orders of 
magnitude in radius. 

Toroidal field instabilities have been suggested to gen- 
erate shocks and high-en ergy emission in quasar jet sys- 
tems JSikora et al J 12005'). For GRB systems, if magnetic 
dissipation continued beyond the 7-ray photosphere at r ~ 
10 7 — 10 8 r- g , then those shoc ks could be directly respon- 
sible for the 7-ray emission (iLvutikov fc Blackmanl 1200 ll : 
iLvutikov fc Blandfordl [2003 ) . This single magnetic dissipa- 
tion model could unify GRB and Blazar emission. This pa- 
per does not resolve the emission mechanism, but the results 
can be used to test the plausibility of emission models. 

Some GRB models assume the Poynting energy 
to be converted into radiation far fro m the collapsing 
star by internal dissipation (see, e.g., |Tliomrjson| Il994j; 
|MeszarofT& Rccs 1997; SDruit^D£d|me ta &^Drenkhahn]|20M|: 
prenkhahnl 120021: iDrenkhahn fc Spruit) 120021: ISikora et*al] 
2003; ILvutikov. Pariev. fc Blandfordl 12003)) such as driven 
by magnetic reconnection. Kink instabilities have been 
invoked to allow reconnection-driven magnetic energy to 
be converted into thermal energy to accel erate the flow 
jDrenkhahrl2002l : IDrenkhahn fc Spruffll2002t) . although the 
reconnection rate assumed in those papers i s not sup- 
porte d by the GEM project on reconnection JShav et alJ 
2001). Some simplified acceleration models assume that kink 
instabilities drive reconn ection that leads to acceleration 
jGiannios fc Spruit) ^006). However, despite the simulated 
jet being marginally stable to current-driven instabilities, 
significant internal dissipation occurs already by 10 3 r 9 . Re- 
connection is not likely nor required. The shocks driven by 
these moderate instabilities are an alternative to the recon- 
nection paradigm. 

We find that after these toroidal field instabilities equal- 
ize the magnetic and thermal energy at r ~ 10 3 r 9 that the 
jet becomes conical for the angular regions with 6j > 5°. 
Consequently, this part of the jet ceases to magnetically ac- 
celerate beyond this radius. In contrast, the core of the jet 
within 6j ~ 5° remains paraboloidal along the entire jet out 
to r ~ 10 4 r 9 . This core may continue to magnetically ac- 
celerate beyond this radius. This suggests that toroidal field 
instabilities play a crucial role in determining the opening 
angle and terminal Lorentz factor of jets. 

Theoretical studies of ideal MHD jets focused 
on cold jets and the conversion of P oynting e nerg 
directly to 
Eichlerl IT99; 



ki netic energy ( s ee, e . 
iBeeel man fc Li 1199- " 



[Li^^^ 



Tomimatsu 



Daiene fc Prenkhahnl 12001 iBeskin fc Nokhrinal I2005T) . In- 

deed, much of the work has focused on self-similar models, 
of which the only self-consis tent solution found is suggested 
to be R-self-similar models <IVlahakisll2004l : iFendt fc Ouvedl 
2004). Unfortunately, such models result in only cylindrical 



© 2006 RAS, MNRAS 000.IHI231 



GRMHD Simulations of Poynting- dominated Jets 19 



asymptotic solutions, which is apparently not what is ob- 
served in AGN jets nor present in the simulations discussed 
in this paper. The previous section showed that some aspects 
of the jet are nearly self-similar and that there is some clas- 
sic ideal-MHD acceleration occurring in this region. How- 
ever, cold ideal-MHD jet models have no way of addressing 
the region where magnetic energy is converted to thermal 
energy in shocks. This region also involves relatively rapid 
variations in the flow, so stationary jet models would have 
difficulty modelling this region. 

6.2 Simulations as Applied to AGN and X-ray 
Binaries 

Quantitatively all of these results can be rescaled by density 
in the jet in order to approximately apply to any GRB model 
and to AGN and x-ray binary systems. This is because the 
jet has only reached about T ~ 10, which is not too small or 
large a Lorentz factor, and Too oc l/pojet- Some astrophysi- 
cal implications of the jet simulations and results as applied 
to AGN are described in the next several paragraphs, while 
a detailed discussion of the results as applied to GRBs will 
be considered in a separate paper. 

The angular structure of the simulated jet may help ex- 
plain many aspects of observed jets. The fast spine can lead 
to strong Compton scattering of photons to higher energy, 
while the outer angular part of the jet can lead to radio 
emission. This process may be required to explain high- and 
low- ener gy emission from TeV BL Lac objects and radio 
galaxies l|Ghisellini et alJl2005l ) and may help explain high- 
energy observations of blazars ( Blandford fc Levinsorlll99,4 
iGhisellini fc Madaul Il996l : IChiaberee et al.ll2000l) . The full 
opening angle of ~ 10° of the simulated jet cor e also agrees 
with the observations of the far-field jet in M87 l|.Tunor et alJ 
119991) . 

The radial structure of the jet and the presence of 
shocks beyond the Alfven surface suggest that synchrotron 
emission should be observed at hundreds of gravitational 
radii from the black hole rather than near the black hole. 
This suggests the shock zone (or "blazar zone" for blazars) 
should be quite extended b etween r ~ 10 2 r Q and up to 
about r ~ 10 4 r 9 (see, e.g.. ISikora et alJ 120051. It is also 
often assumed that if the jet is highly collimated that it is 
also highly relativistic near the black hole, which would sug- 
gest Comptoniza tion of disk photons sho uld produce clear 
spectral features (|Sikora fc Madeiskil200of) . However, the jet 
may rather accelerate slowly but collimate quickly, which is 
what we find. There is an early transfer of Poynting flux 
to kinetic energy flux leading up to about r ~ 5 — 10 by 
about r ~ 10 3 r 9 . This is consistent with the lac k of observed 
Comp tonization features in blazars (see, e.g., ISikora et alJ 
2005). At large distances the field becomes toroidally domi- 
nated, but this does not necess arily contradict observations 
of pa rallel fields in FRI sources (|Blandfordl2000l : ISautv et alJ 

Prior work suggested that the BZ power is insufficient t o 
account for Blazar emission l|Maraschi fc Tavecchiolfeooffi . 
where they assumed that T ~ 10 and 8j ~ 15°. However, 
we find that the structure of the Poynting-dominated jet is 
nontrivial. The region with T ~ 10 is narrower with 9j ~ 5° 
and jet emission may be dominated by shock accelerated 
electrons with thermal r' MA -* ~ V > 100 with an extended 



high energy tail. This lowers the necessary energy budget of 
the jet enough to be consistent with BZ power driving the 
jet. 

The disk wind may play a crucial role independent of 
the jet. For example, relatively thin disks or slowly rotating 
black holes wou ld produce disk winds that could appear as 
"aborted jets" iGhisellini et alJ 12004) The classical AGN 
unification models JUrrv fc Padovanilll995l) invoke a dom- 
inant role for the molecular torus and broad-line emitting 
clouds, while the broad disk wind composed of magnetized 
blob-like regions may significantly contribute to mo difica- 
tions and in understan ding the origin of the clouds jElvisI 
l2000UElvis et alJl2004) . 

Erroneous conclusions could be drawn regarding the 
jet composition since there is actually a jet and disk 
wind, which could each have a separate composition. That 
there can be two separate relativistic jet components 
makes it diffic ult to draw c l ear c o nclusions regarding the 
composition ([G uilbert et alj[l98^: [Celotti fc Fabianl Il993l : 
iLevinson fc Blandfordlll996t ISikora fc Madeiski 120001 . En- 
trapment, which could occur at large distances when the 
ideal MHD approximation breaks down, also causes difficul- 
ties in isolating the "proper" jet component's composition. 
In cases where only electron-positron pair jets can be sup- 
ported, then this would indicate that disk winds are ruled 
out for such sources, such as r ecently suggested for FRII 
sources l|Kino fc Takaharalf2004 l . 

6.3 Disk Winds vs. Poynting Jets in M87 

M87 is associated with the most widely studied AGN jet, so 
here we discuss implica tions of our results for M87 sp ecifi- 
cally. I.Tunor et alJ l|l99flh and lBiretta et alJ jl99fltl2002n sug- 
gested that M87 slowly collimates from a full opening angle 
of about 60° near the black hole to 10° at large distances. 
However, some of their assumptions are too restrictive. First, 
they assumed t he jet is always co nical, which is apparent 
from figure 1 in ljunor et alJ il999ll . If the jet is not conical 
this can overestimate the opening angle close to the core (i.e. 
perhaps 35° is reasonable all the way into the core). Second, 
their beam size was relatively large so that factors of 2 error 
in the collimation angle are possible. Third, they assumed 
the radio core is the location of the black hole, while the 
inner-jet may be cold and only produce synchrotron emis- 
sion once shocks develop. Finally, the assumption that the 
jet is a single, s lowly collimating jet remains t he prevailing 
view (see, e.g., iTsinganos fc Bogovalovl 120051 iGracia et alJ 
2005). GRMHD numerical simulations of jet and disk wind 
formation suggest that a model with a single jet is too sim- 
ple of an interpretation of the observations. If there is a 
highly collimated relativistic Poynting jet surrounded by a 
weakly collimated disk wind, then this would also fit their 
observations without invoking slow collimation. 

A form of the idea that winds are coincident 
with, and collimate, jets h as also been proposed by 
ITsinganos fc BogovalovM2005l) and applied to M87, but they 
considered a model where the wind slowly collimates the jet 
in order to fit observations. Here we suggest that the obser- 
vations may have been misinterpreted due to the presence 
of two components: a well-collimated relativistic cold Poynt- 
ing jet and a mildly relativistic disk wind. We suggest the 
broader emission component could be due to the disk wind. 
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More recent maps of the M87 jet formation region 
show ed no "jet formation" structure (iKrichbaum et all 
I2004T) . Thus, the structures seen previously may be tran- 
sient features, such as associated with turbulent accretion 
disk producing a dynamic disk wind or associated with time- 
dependent shocks produced in the jet far from the black hole. 

Measurements of the apparent jet speed in M87 reveal 
typically T ~ 1.8 near the core while F ~ 6 at larger radii. 
However, so me core regions are associated with F ~ 6 that 
rapidly fade feiretta et al.lfl998l) . This is consistent with a 
two-component outflow where the cold fast moving core of 
the jet is only observed if it interacts with the surround- 
ing medium, the slower disk wind, or it undergoes internal 
shocks. 

6.4 Comparisons with Other Numerical Work 

Numerical models of jets continue to be studied using the 
hydro d ynamic ap p roxim ati on (iMa^F^^gri^y^c^sJ^vl 
19991: lAlovetal.1 Eoool. 120021: IScheck et all 120021: 

Zhang. Wooslev fc MacFadvenl Bool IZhang W. et alJ 
20041: lAlov et al.ll2005h . despite the presence of dynamically 
important magnetic fields in the accretion disks that 
have been shown to lead to strongly magnetized jets 
llDe Villiers. Hawlev. fc KroliklEool : iMcKinnev fc Gammiel 
l2004l:lDe Villiers et alJl2005aT) . These hydrodynamic studies 
have suggested that hydrodynamic mixing within the 
jet may lead to observational signatures in the emission 
variability. Other models suggested that such hydrodynamic 
mixing leads to an inability to clearly distin guish the jet 
comp osition as either leptonic or baryonic ijScheck et alJ 
2002). We suggest magnetic instabilities due to pinch or 
kink modes dominate the variability rather than hydro- 
dynamic instabilities, which are kno wn to be qu e nched 
by magnetic con f inement (see, e.g., |Clarkee^^J_^8d 
i Lind et all Il989t IAddI fc Camenzindl Il992l : iRosenet^dl 
1999). The magnetic confinement severely limits the mixing 
between the jet and the environment. 

Hydrodynamic studies have suggested that shocks 
driv en into the surrounding environment collimate the 
iet jAlov et al]|200rilzhang. Wooslev. fc MacFadvenl l2003t 
IZhang W^t^lT2004T) . We find that the jet structure is neg- 
ligibly broader for models with a negligible surrounding en- 
vironment. This is because the jet is magnetically confined 
and collimated at large radii. As applied to GRBs, our re- 
sults suggest that the presence of an extended stellar enve- 
lope plays no role in the jet dynamics except to provide the 
disk with material. However, the jet plays a significant role 
in modifying the stellar envelope by sending lateral shocks 
through the star. 

Similar studies have also suggested that variations 
in the Lorentz f actor are due to hydrodynamic effects 
iAlov et alJl200 j . 120051) . but we find that variations in the 
Lorentz factor are due to magnetic instabilities. Unlike their 
models, our models generate a relatively simple jet structure 
that can be fit by a Gaussian or top-hat + exponential wings. 
Once the shock heating generates an equilibrium magnetic 
fireball, the internal thermal support and toroidal magnetic 
confinement keep the jet stable and conical out to large dis- 
tances. For an optically thick flow like in GRB jets, the flow 
should convert the thermal energy to kinetic energy in adi- 
abatic expansion and lead to F ~ 10 3 . 



For hydrodynamic models of GRBs, en ergy is in- 
jected to model the annihilation of neutr inos ||AloY_et_al 
1 2000 1 : IZhang. Wooslev. fc MacFadvenll2003l : IZhang W. et al' 
2004). They assumed a highly relativistic jet is formed early 
near the black hole and they injected the matter with a large 
enthalpy per baryon of about u/poc 2 ~ 150. They tuned 
the inj ected energy to h ave their results agree with observa- 
tions dFrail et aljboOll) . Other models with a disk and jet 
find that u/poc 2 ~ 10, which is sug gested to imply that 
Fqc ~ 10, insufficient for GRBs jMacFadven fc Wooslev! 
1999), see their figure 28. In other work we suggested that 
the injected energy per baryon is only u/poc 2 < 20 unless 
super-efficient neutrino mechanisms ar e invoked dMcKinnevI 
2005bl: iRamirez-Ruiz fc Socra tes 2005). We also found that 
the energy released is larger than observed, suggesting an 
fairly inefficient generation of 7-rays. We find an energeti- 
cally dominant, lower F jet component that may explain x- 
ray flashes. In our case the baryon-contamination problem is 
avoided by magnetic confinement of the jet against baryons 
from the disk. We find F x ~ 100 — 1000, which is much 
larger than they found and this may avoid the compactness 
problem. The difference between their and our model is the 
presence of a magnetic field and a rotating black hole, which 
together drive the BZ-effect and a stronger evacuation of the 
polar jet region. 

MHD simulations of jet propagation have been per- 
formed that correspond to weakly magnetized AGN jets , 
where b 2 /{2u) < 3.3 and b 2 /p < 0.2 jLeismann et al.l2005ll . 
Our study suggests that this region may only be important 
for AGN Poynting-dominated jets at scales r 3> 10 4 r 9 once 
the jet becomes matter-dominated or may apply to the disk 

outflow. 

Nonrelativistic MHD dProga et al] l2003t iKato et all 

| 2004l) and prior G RMHD simulations jMizuno et al]l2004l : 
lKomissaro\ll2005ll of jets showed a slow jet with v ~ 0.2 — 
0.3c, insufficient to explain most jets from GRBs and AGN. 
iProga et al] (120031) suggested that their Poynting-dominated 
jet has Too < 10, but they were unable to follow the rela- 
tivistic jet using the nonrelativistic equations of motion. The 
previous GRMHD simulations created slow jets because of 
the "floor model" in the polar regions, the short time of in- 
tegration, and the limited outer radius of the computational 
box . 

iDe Villiers et al ] j2005bTl used a nonconservative nu- 
merical method to evolve a fully relativistic, black hole mass- 
invariant model, and showed a jet with hot blobs moving 
with r ~ 50. A primary result in agreement with the results 
here is that a patchy or pulsed "magnetic fireball" is pro- 
duced. This suggests that the development of a "magnetic 
fireball" is not an artifact of the numerical implementation 
but is a result of shock heating. We also agree in finding 
that the core of the jet is hot and fast and is surrounded by 
a cold slow flow. One difference is that they say they seem 
to have found that the flow is cylindrically collimated by 
r > 300r 9 , while we find nearly power-law collimation un- 
til r ~ 10 3 r 9 and an oscillatory conical asymptote beyond. 
They suggested that temporal variability is due to injection 
events near the black hole, while we suggest it is due to pinch 
(or perhaps kink) instabilities at r > 10 2 r 9 . 

They found a larger value of F at smaller radius than we 
find. Indeed, one should wonder why their black hole mass- 
invariant model applies to only GRBs and not AGN or x-ray 
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binaries. This is a result of their much lower "floor" density 
of p Q ~ 10~ 12 po, disk, which is far too low to be consistent 
with the self-consistent pair-l oading l|Pop ham et al]ll9 99h 
or neutron diffusion loading llLevinson & Eichlcr 2003) of 
the jet. Also, because they used a constant density floor, 
the jet region at large radius must be additionally loaded 
with an arbitrary amount of rest-mass. The typical "floor" 
model adds rest-mass into the comoving frame, but this ar- 
tificially loads the jet with extra mass moving at high V. 
As applied to GRBs, we suggest that rather than the ac- 
celeration occurring within r < 700r 9 of their simulation, 
that acceleration occurs much farther away before the emit- 
ting region at r < 10 9 r 9 . In their F ~ 50 knots the gas has 
10 3 ^5 u/poc 2 < 10 6 . This implies that their actual terminal 
Lorentz factor is on the order of 10 5 < T x < 10 7 , which 
would imply that the external shocks occur before internal 
shocks could occur. 

6.5 Limitations 

As has been pointed out by iKomissarovl J2005T) , noncon- 
servative GRMHD schemes often overestimate or underesti- 
mate the amount of thermal energy produced. All numerical 
models suffer from some numerical error. Shock-conversion 
of magnetic energy to thermal energy is modelled by the 
conservative scheme in HARM in the perfect magnetic fluid 
approximation with total energy conserved exactly. How- 
ever, our numerical model may still have overestimated the 
amount of magnetic dissipation in shocks, and a more accu- 
rate calculation may more slowly convert magnetic energy 
to thermal energy. However, we expect that toroidal field 
instabilities drive efficient magnetic dissipation as shown in 
the numerical results presented here. It is also encouraging 
that s imilar result are fo und in 3D nonrelativistic simula- 
tions JOuved et alJl2003h and maybe in the 2D relativistic 
simulations of iDeViffiers et alJ i2005bTl (they did not say 
why their jet generates hot blobs) . 

Co nfidence in the numer ical results is obtained by code 
testing ( Gammic et all2 003a). comparisons w ith related an- 
alytical models (iMcKinnev fc Gammiell2004D . and conver- 
gence testing the specific model being studied. For the fidu- 
cial model simulated in this paper, we used both higher and 
lower order spatial and temporal reconstructions and used a 
resolution smaller by a factor of two. The results described 
in the paper are independent of these numerical changes. 

Unlike GRB neutrino-dominated disks, AGN and X- 
ray binaries may have accretion disks with a wide variety of 
H/R near the black hole, which could lead to a wide to nar- 
row opening angle of the jet and disk wind. This broadness 
of the jet and wind may affect the classical view ing- angle- 
dependent unification models JSautv et alj l2*002). Only a 
self-consistent radiative GRMHD calculation can determine 
the disk thickness. The quantitative conclusions in this pa- 
per regarding the collimation angle assumed H/R ~ 0.2 near 
the b lack hole, while H/R ~ 0.9 fADAF-like. lNaravan &y] 
Il995l) is perhaps more appropriate for such systems. The 
sensitivity of these results to H/R is left for future work. 

The simulations used a field geometry that was initially 
moderately organized, although there was no net field. As 
discussed in the introduction and in the floor model sec- 
tion, moderate changes in the initial geome try do not af- 
fect the results dMcKinnev fc Gammiell2004h . However, ac- 



cretion of a highly irregular (tangled) field would be quite 
different. An organized field may not develop around the 
black hole, and so disk material would mass-load the jet. 
In thi s case, the jet just becomes an extension of the disk 
wind dMcKinnev fc Gammiell2004h . However, the existence 
of a mostly uniform field threading the disk arises naturally 
during core-collapse supernovae and NS-BH collision debris 
disks. In AGN and stellar wind- capture x-ray binary sys- 
tems, the accreted field may be uniform over lon g time scales 
jNaravan et alJl2003l: ISpruit fc Uzdenskvl2005l) . Roche-lobe 
overflow x-ray binaries, however, may accrete quite irregu- 
lar field. The field geometry that arrives at the black hole, 
after travelling from the source of material (molecular torus, 
star(s), etc.) to the black hole horizon, should depend sen- 
sitively on the reconnection physics. 

In order to evolve for a longer time than simulated in 
this paper, other physics must be included. For long-term 
evolution of a GRB model, one must include disk neutrino 
cooling, photodisintegration of nuclei, and a realistic equa- 
tion of state. If one wishes to track nuclear species evolution, 
a nuclear burning reactions network is required. For the neu- 
trino optically thick region of the disk, radiative transport 
should be included. The self-gravity of the star should be in- 
cluded to evolve the core-collapse. This includes a numerical 
relativity study of the collapse of a rotatin g magnetized mas- 
sive star into a black hole (see review bv lStergiouIasll2003l . 
§4.3). The jet should be followed through the entire star 
and beyond penetration of the stellar sur face l|Aio^_et_al 
| 2000tlZhang. Wooslev. fc MacFadvenll2003l : IZhang W. et~ 

As applied to all black hole accretion systems, some 
other limitations of the numerical models presented include 
the assumption of axisymmetry, ideal MHD, and a nonra- 
diative gas. 

The assumption of axisymmetry is not crucial for the 
basic structure of the inner je t region since our earlier results 
jMcKinnev fc G ammie||2004j) a re in quantitative agre ement 
with 3D results JPe Villiers. Hawlev! fc Krolikl 120031) . The 
primary observed limit ation of axisym metry appears to be 
the decay of turbulence dCowling|l934l) . which we attempted 
to avoid by requiring a resolution that gives quasi-steady 
turbulence for much of the simulation. Also, the jet at large 
distances has already formed by the time turbulence decays, 
and by that time the jet at large radius is not in causal 
contact with the disk. 

When the toroidal field dominates the poloidal field, 
eventually m = 1 kink instabiliti es and higher modes 
may appear in 3D models (see, e.g ., INaka mura et al.ll200ll : 
lOuved et a!1l2003l : INakamura fc Meierll2004l) . Thus our 2D 
axisymmetric models may have underestimated the amount 
of oscillation in the flow and the conversion of Poynting 
flux to enthalpy flux. On the contrary, rotating Poynting- 
dominated jets may be m arginally stable to kink instabilities 
llTomimatsu et al] l2001). as suggested by our results. 

The decay of turbulence due to the limitation of ax- 
isymmetry may also limit the efficiency of the Blandford- 
Znajek process. The magnetic arrested disk (MAD) model 
suggests that any accretion flow that accumulates a large 
amount of magnetic flux near the black hole eventually halts 
the accretion flow and bu i lds a black hole mag netosphere 
ilgumenshchev et al]l2003l ; iNaravan et al . 2003). This im- 
plies that the efficiency of extracting energy could be higher. 
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We have neglected high-energy particle and radiative 
processes in the numerical simulations. The floor model of 
the minimum allowed rest-mass density and internal energy 
density is ad hoc and chosen so that a fast jet emerges. Fu- 
ture work should include a model of pair creation and pair 
annihilation in order to simulate self-consistent mass-loading 
of the Poynting-dominated jet. A self-consistent model of the 
mass-loading can at least determine the density that should 
be present near the black hole around the poles. However, 
accurate determination of the Lorentz factor requires at least 
radiative transfer and Comptonization to model the radia- 
tive reaction within the jet that can slow or speed the jet. 

For AGN and x-ray binaries, the radiatively inefficient 
disk approximation, which assumes electrons couple weakly 
to ions, may not hold. If the electrons and ions eventually 
couple near the black hole, then the disk might collapse into 
an unstable m agnetically dominated accretion disk (MDAF) 
jMeieill2005l) . Like the MAD model, this might drastically 
alter the results here, although it is uncertain whether jets 
are actually produced under the conditions specified by the 
MDAF model. 

The single-fluid, ideal MHD approximation breaks 
down under various conditions, such as during the quies- 
cent output of AGN and x-ray black hole binaries, where a 
two-temperature plasma may form n ear the black hole a s 
ions and electrons decouple (see, e.g.. iNaravan yI IiQQS). 
In some AGN and x-ray binary systems, the gas may no 
longer act like a fluid and nonthermal Fermi acceleration 
can produce a jet in a shear acceleration reg ion between 
the disk and corona dSubramanian et al.lll999T) . Resistivity 
plays a role in current sheets where reconnection events may 
genera te flares as on the sun, such as possibly observed in 
Sgr A jGenzel et al.ll2003l) . Finally, radiative effects may in- 
troduce d ynamically important instabilities in th e accretion 
disk fe.g. lGammilll99dlBiaes fc Socratesll2003D . 



7 SUMMARY 

A GRMHD code, HARM, was used to evolve black hole 
accretion disk models. Our work extends the results of pre- 
vious GRMHD numerical models by studying the Poynting- 
dominated jet in more detail and studying the jet out to 
r ~ 10 4 r 9 till t ~ 10%. 

Basic results and conclusions include: 

(i) Poynting jet reaches r ~ 10 following nearly 
paraboloidal field lines and nearly F oc r 1//2 . Beyond r ~ 
10 3 r 9 , the jet has a conical outer angular part and a nearly 
paraboloidal core. 

(ii) Peak energy flux within the Poynting jet goes from 
6j ~ 60° near the black hole to 0j ~ 5° at large distances. 

(iii) Energy structure of the Poynting jet is Gaussian with 
a half- width of 9 « 8° . 

(iv) Poynting jet has a core with flat energy flux within 

(v) Extended slower Poynting jet component with 6j ~ 
27°. 

(vi) Poynting jet variability is dominated by toroidal field 
instabilities. 

(vii) Poynting jet is marginally kink stable. 

(viii) Poynting flux is shock-converted into enthalpy flux 
beyond Alfven surface in an extended "shock zone." 

(ix) Maximum terminal Lorentz factor is Too < 10 3 . 

(x) Poynting flux and enthalpy flux come into equiparti- 
tion by r ~ 10 3 r 9 . 

(xi) For radiation-dominated conditions, optically thick 
flows can tap thermal energy in adiabatic acceleration, while 
optically thin flows can lose thermal energy as radiation. 

(xii) Disk wind has collimated edge (near the Poynting 
jet) with T < 1.5. 

(xiii) Full disk wind is broad even near the black hole, 
and this may account for significant spatially-broad emission 
from near AGN cores, such as in M87. 



6.6 Future Work 

As applied to the collapsar model and other GRB models, 
future calculations will include a realistic equation of state; 
photodisintegration of nuclei; gene ral relativistic ray-tracin g 
neutrino transport (similar to, e.g., Brodcrick & Loeb 2005); 
neutrino cooling (similar to, e.g.. lKcliri^t^lTl2^05l) ~ gen- 
eral relativistic accounting for the neutrino annihilation; 
neutrino Comptonization; simplified photon transport ; and 
photon Comptonization. The numerical calculations will 
also be performed in 3D. 

As applied to AGN and x-ray binaries, future calcula- 
tions will in clude general relativistic photon transport (sim- 
ilar to, e.g., iBroderick fc Loebll2005l) : photon cooling; gen- 
eral relativistic accounting for the photon annihilation in the 
Poynting-dominated jet; photon Comptonization; and ther- 
mal and nonthermal synchrotron emission. The numerical 
calculations will also be performed in 3D. 

The same type of calculation can be performed for sys- 
tems harboring neutron s tars, and this is a n atural extension 
of recent force- free work <lMcKinnevll2006al lbl) . 
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